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ABSTRACT 

The  goal  of  this  thesis  was  to  develop  and  validate  a  set  of 
microcomputer  programs  using,  in  part,  a  previously  written 
finite-element  algorithm  to  analyze  the  perturbation  of  an  incident 
electromagnetic  field  as  it  penetrates  a  missile  radome.  An 
interactive  program  was  developed  to  design  the  radome  structure 
using  aerodynamic  shaping  functions  and  provide  structure  and 
material  files  as  input  to  the  finite-element  code.  A  second  program 
was  developed  to  use  the  spherical  harmonic  expansion  coefficients 
provided  by  the  finite-element  code  to  assemble  the  electromagnetic 
fields  within  the  radome  core  and  to  display  the  fields  that  appear 
across  a  planar  antenna  with  three-dimensional  graphics  for  any 
orientation  of  the  antenna.  Algorithms  were  included  which  compare 
the  computed  field  components  to  the  theoretical  incident  plane 
wave  for  each  stipulated  angle  of  incidence,  in  order  to  determine 
the  perturbation  due  to  the  presence  of  the  radome.  Validation  of 
the  computational  method  was  attempted  by  analyzing  the 
perturbation  indicated  for  an  ideal  radome  with  relative  permittivity 
of  unity.  The  attempted  validation  showed  phase  errors  in  the 
computed  fields  which  were  minimal  for  axial  incidence,  but  became 
significant  for  highly   canted   incidence. 
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I.  INTRODUCTION 

Radome  design  is  a  composite  discipline  evolved  from 
aeronautical,  electrical,  mechanical  and  materials  engineering.  The 
word  radome  was  coined  during  World  War  II  as  a  composite  of  two 
words,  radar  and  dome  [Ref.  l:p.  ix].  The  term  radome  originally 
described  the  family  of  dome-shaped  structures  designed  and  built 
to  house  and  protect  the  radar  antennas  on  aircraft,  but  now  the 
term  is  applied  to  any  structure  that  protects  a  radio  frequency 
antenna  system.  The  structure  is  designed  to  provide  maximum 
protection  to  the  antenna  with  minimum  interference  to  the 
antenna's    electrical    performance. 

Radomes  on  missiles  are  presently  designed  to  meet  many 
diverse  specifications.  These  specification  requirements  can  be 
divided  into  two  categories,  electrical  and  aero-mechanical.  Electrical 
design  problems  include  those  that  have  to  do  with  the  radome's 
ability  to  act  as  a  window,  ideally  transparent,  to  the  electromagnetic 
energy  that  must  pass  through  the  radome  to  be  used  by  the 
missile's  antenna  system.  This  energy  is  of  a  specific  frequency 
range  and  intensity  that  is  within  the  operating  parameters  of  the 
antenna  system.  The  magnitude  and  form  of  the  electromagnetic 
energy  will  be  dependent  on  range  from  the  source,  polarization,  and 
incident  angle.  Aero-mechanical  design  problems  address  not  only 
the  design   of  the   structure   of  the   radome,   but   also   the   ability   of  the 
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radome  to  withstand  environmental  conditions  and  battle  damage. 
Environmental  conditions  vary  with  scenario,  and  include:  normal 
pressure  force  and  viscosity,  wind  loads,  inertial  loads,  temperature 
variations,  moisture  resistance,  resistance  to  rain  erosion,  and  other 
diverse  factors.  Design  of  a  radome  for  use  in  battle  may  require 
that  the  aero-mechanical  design  consider  resistance  to  nuclear 
radiation    and    laser   weapons. 

As  in  any  engineering  design,  there  are  tradeoffs  between  the 
various  requirements,  such  that  each  interrelated  component  design 
factor  must  be  weighted  and  considered  in  the  overall  system  design. 
Numerical  models  of  environmental  conditions  assist  in  the  analysis 
of  the  tradeoffs  for  the  aero-mechanical  design.  Numerical  models  to 
aid  in  the  electrical  design  analysis  of  a  radome  to  protect  missile 
antenna  systems  which  operate  at  high  frequencies  are  available, 
such  as  the  ray-tracing  algorithms  described  by  Hirsch  and  Grove  in 
Ref.  2.  When  a  missile  antenna  system  requires  operation  in  the 
resonance  region  of  frequencies  for  the  radome  parameters,  a 
numerical  model  is  required  that  can  solve  the  radome  interaction 
problem. 

The  subject  of  this  thesis  is  the  development  of  a  generalized, 
interactive  set  of  programs  to  be  used  in  the  analysis  of  the 
interaction  between  an  actual  or  proposed  radome  structure  and 
incident  electromagnetic  fields  at  the  lower  frequencies  of  interest, 
where  strong  resonance  characteristics  exist.  In  order  to  analyze  a 
specified  radome,  a  numerical  model  is  constructed  using  the 
interactive    program,    RADSHP.        This    radome    generation    program 
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provides  geometric  shaping  functions  to  allow  the  user  to  construct  a 
multilayer  radome  model  of  specified  materials.  Chapter  II  provides 
a  discussion  of  radome  design  methodology,  the  aerodynamic  shaping 
functions  available  in  this  package,  radome  wall  constructions,  and 
material    specification. 

This  thesis  work  is  an  application  of  the  programs  developed 
for  use  on  a  personal  computer  by  Morgan  [Ref.  3]  and  Connolly 
[Ref.  4]  to  calculate  electromagnetic  scattering  by  inhomogeneous 
dielectric  and/or  magnetic  bodies  of  revolution.  As  discussed  in 
Chapter  II,  radomes  are  manufactured  using  dielectric  materials  in 
order  to  provide  a  radio  frequency  window  to  the  antenna  housed  in 
the  radome's  interior.  An  axisymmetric  missile  radome  structure 
meets  the  physical  structure  requirement  for  rotational  symmetry. 
The  numerical  analysis  of  the  interaction  of  the  specified  radome 
structure  and  the  incident  electromagnetic  energy  is  performed  by 
the  code  developed  by  Morgan,  EMCAD.  Modifications  to  this  code 
resulted  in  the  program  EMRAD,  which  provides  data  file  generation 
to  store  the  modal  expansion  coefficients  which  are  used  by  the 
program  CORFLD  to  determine  the  electromagnetic  fields  within  the 
radome  interior.  In  Chapter  III,  the  theory  required  for  the 
formulation  of  the  computed  interior  fields  of  the  radome  by  the 
program  CORFLD  is  discussed. 

The  programs  developed  by  Morgan  use  an  optimized 
variational  finite-element  algorithm  in  conjunction  with  a 
tri-regional  unimoment  method  to  provide  scattering  solutions  for 
each    of  multiple   incident   fields   impinging   on    the    specified    structure. 
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The  programs  represent  a  significant  improvement  in  computational 
speed  and  versatility  over  computer  codes  which  find  the  solution  to 
the  electromagnetic  resonance  behavior  by  use  of  surface  current 
integral  equation  formulations.  An  overview  of  the  finite  element 
scattering  algorithm,  along  with  the  concepts  embodied  in  the 
coupled  azimuthal  potentials  and  the  unimoment  method  are 
included  in  Appendix  A.  Further  information  is  available  in 
references  3,  4,  5,  6,  and  7. 

Chapter  IV  discusses  the  flow  of  data  between  the  programs 
which  have  been  developed.  The  general  use  of  the  RADSHP 
program  to  construct  the  radome  data  files  is  discussed,  as  well  as 
the  data  files  generated  and  their  uses.  The  flow  of  information  into 
and  out  of  the  radome  analysis  program  EMRAD  is  analyzed.  The 
execution  and  operation  of  CORFLD  is  discussed. 

In  Chapter  V,  the  validation  of  the  interior  field  computation  is 
discussed,  wherein  the  numerically  generated  fields  are  compared  to 
the  incident  fields  for  a  radome  constructed  of  air.  The  results 
produced  for  radomes  constructed  of  materials  presently  used  in  the 
radome    industry    are   presented    and    discussed. 

The  conclusions  in  Chapter  VI  will  summarize  the  work 
previously  presented,  as  well  as  provide  recommendations  for  future 
modifications   to  the  interactive  package  of  programs   being   used. 

Additional  Appendices  are  included  to  provide  copies  of  the 
source  code  for  the  computer  programs,  RADSHP  and  CORFLD,  which 
were  written  as  part  of  this  thesis  effort.  Source  listings  for 
additional    programs,    such    as    EMCADIN,    EMESH,    and    EMRAD    are 
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available   upon  request  from  Prof.   Michael   A.   Morgan.     His   address  is 
given  in  Appendix  D,   Software  Sources. 


II.  RADOME  DESIGN 

A.  RADOME  DESIGN  METHODOLOGY 

The  configuration  of  a  radome  for  high-speed  missiles  must 
combine  high  aerodynamic  performance  with  a  highly  accurate 
transmission  of  electromagnetic  energy  through  the  radome  to  the 
missile  seeker.  As  described  in  Ref.  2,  the  methodology  for 
establishing  radome  performance  requirements  is  an  iterative 
process.  Figure  2.1  shows  the  process  as  a  set  of  design  procedure 
loops  that  begin  with  a  statement  of  desired  performance,  and 
include  analytical  studies,  hardware  development,  and  hardware 
evaluation.  These  result  in  a  set  of  optimized  radome  requirements, 
which  is  the  blueprint  from  which  a  radome  can  be  constructed. 
Computer  software  tools  help  to  maximize  the  benefit  of  the  work 
done  in  analytical  studies,  thereby  minimizing  the  number  of 
iterations  required  for  the  hardware  development  and  evaluation. 
The  programs  described  in  this  thesis  are  designed  to  analyze  the 
interaction  between  the  aerodynamic  considerations  in  radome 
design   and   the  electrical  performance  of  the  radome. 

The  aerodynamic  considerations  in  the  design  of  a  radome 
structure  are  three-fold:  shape-generated  drag,  thermal  stress  from 
aerodynamic       heating,       and       structural       performance.       These     three 
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Figure   2.1      Methodology   for   establishing   radome   requirements 


parameters  are  not  independent  of  each  other,  and  the  electrical 
performance  of  the  radome  is  a  strong  function  of  these  aerodynamic 
factors.  The  outer  shape  of  the  radome  has  a  large  impact  on  how 
the  missile  will  fly  through  the  air.  It  is  this  outer  conformation  that 
will,  along  with  speed  of  travel  and  angle  of  attack,  determine  the 
forces  that  act  on  the  radome  structure  during  flight.  The  forces 
generated  during  flight  determine  the  design  requirements  on  the 
radome  structure  for  material  strength  and  thermal  stress. 
Fortunately,  the  shapes  that  result  in  the  best  drag  performance  are 
those  that  can  also  be  used  to  develop  radomes  that  match  the 
structural  requirements.  Unfortunately,  the  shapes  that  best  meet 
the  drag  requirements  are  not  necessarily  the  best  for  maximizing 
the   electrical   performance.  The   programs   developed   in   this   thesis 

include  the  computer  aided  design  (CAD)  functions  to  analyze  several 
shaping  functions  that  aerodynamicists  have  developed  to  minimize 
the  aerodynamic  drag  forces  with  respect  to  radome  structure 
parameters. 

B.  RADOME  SHAPING  FUNCTIONS 

An  ogive  shaping  function  is  one  of  the  most  common  used  in 
radome  designs.  The  ogive  shape  is  popular  because  it  is  easily 
fabricated  and  has  relatively  acceptable  aerodynamic  performance  at 
high  speeds.  It  is  defined  [Ref.  l:p.  48]  as  a  segment  of  the  arc  of  a 
circle,    whose    radius,    Fj,    is    greater   than    the   radius    of   the    radome 

base.      A    pointed   ogive   is    generated    when    the   circle's   radius,   Fj,  is 
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large   enough   so   that   the   arc   reaches   the   radome's   axis   of  symmetry. 
The  ogive   shaping   function   is  described  by 


=     1/f12-((H-z)-F2)2    -  F3   , 


(2.1) 


where  F2  and  F3  determine  the  position  of  the  circular  shaping  arc,  as 
illustrated  in  Fig.  2.2.  The  ogive  shaping  function  degenerates  to  a 
hemisphere     when     the     shaping  factors,     F2    and  F3,   are  defined   as 

F2  =  Fj    and  F3  =  0. 

A  tangent  ogive  shaping  function  generates  a  radome  shape 
such  that  the  surface  of  the  radome  is  tangent  at  the  base  of  the 
radome,  providing  a  smooth  transition  to  the  missile  body  [Ref.  8].  A 
tangent   ogive   shaping   function   is   determined    such   that 


x   =  4 


JFa-J    ~  Ff 


(2.2) 


where 


F,   = 


*  +  2Mi 

2  R 


(2.3) 


and 


F<    =     F.  -    ^ 


(2.4) 


A     shaping    function     that    provides    maximum    volume-to-drag 
ratio  at  high  flight  speeds  is  the  von  Karman  function  [Ref.  l:p.  48]  , 
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x   = 


rVti 


E   - 


-sin(2Fz) 


(2.5) 


z-axis 


Figure     2.2       Radome  shape  -  nomenclature 
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where 

(H-z) 


Fz   =     cos 


-l 


1   -   2 


(2.6) 


H 

Radomes  described  by  this  von  Karman  shaping  function  have  a 
limited  application  due  to  fabrication  difficulties  related  to  the 
complex    shape. 

In  order  to  minimize  drag  for  a  given  height-to-radius  ratio,  a 
Newtonian  contour  is  used.  This  contour  is  approximated  by  the 
power  series  radome  shaping  function   [Ref.   l:p.  49], 


x  =      R 


'(H-z)^6 


(2.7) 


H 

The   exponent,    F6,  provides  for  shaping.     When  F6  =   1.0,  the  radome 

shape   generated   is   conical.      Aerodynamicists   normally   use   F6  =  0.75 

to   approximate   the   Newtonian   contour. 

A  shaping  function  analyzed  by  Hirsch  and  Grove  [Ref.  2]  uses  a 
parabolic   basis   function   for  the   radome   shape, 

x2   =  F7-z  .  (2.8) 

The  rounded  effect  developed  at  the  tip  of  the  radome  by  the 
parabolic  shape  function  can  be  eliminated  by  the  creation  of  a 
mathematical  "deadzone"  for  points  within  a  region  about  the 
radome's  apex.       The  parabolic  shaping  coefficient,  F7  ,  determines  the 

extent    of    the    "dead    zone"    and    thereby    controls    the    slope    of    the 

parabolic    surface.      The   radome    shape    is    determined    in    terms    of   the 

shaping   coefficient,    then    scaled   to   the   radome   size   parameters.      The 
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unsealed    shaping    function, 

x  =      ^-z    ,  (2.9) 

is  used  to  generate  a  parabolic  radome  shape  with  a  height  of  F7  .  In 
order  to  eliminate  the  parabolic  rounding  of  the  apex,  the  "deadzone" 
is  established  for  any  points  for  the  unsealed  z  value  within   1.0  of  F7. 

Within  this  deadzone,  z  >  (  F7  -  1  ),  the  radome  shape  is  determined 
at  a  line  of  constant  slope,  converging  at  a  height  of  z  =  F7  +  1.0.  This 
results  in  the  marrying  of  a  pointed  apex  onto  the  parabolic  radome 
shape,   with   the  magnitude  of  the   shaping  coefficient,   F7  ,    determining 

the  relative  magnitude  of  the  apex   to  the  total   radome   height. 

This  project  contains  coding  in  RADSHP  to  generate  radomes  of 
multiple  layers  for  the  shaping  functions  described:  the  general 
ogive,  the  tangent  ogive,  the  von  Karman  shape,  the  power  series 
radome,  and  the  parabolic  radome  with  pointed  apex.  RADSHP  can 
generate  a  complex  radome  construction  model  using  different 
shaping  functions  for  each  layer.  RADSHP  generates  a  material  data 
file  appropriate  for  input  to  the  numerical  analysis  program,  EMRAD, 
and  a  structure  data  file  that  can  be  viewed  and  modified  using 
CURVE  DIGITIZER  (Appendix  2)  ,  then  translated  by  EMCADIN  into  an 
input  structure  data  file  for  EMRAD. 

C.  RADOME  WALL  CONSTRUCTION 

If  a   material   existed    which    was   mechanically    strong   enough    to 

sustain    the    aerodynamic    forces    that    the    radome    experiences    during 
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flight,  and  that  same  material  had  the  electromagnetic  properties  of  a 
vacuum,  then  the  radome  could  be  built  to  ideally  meet  all 
requirements.  Such  a  material  does  not  exist,  and,  in  general,  as  the 
requirement  for  material  strength  increases,  the  electromagnetic 
properties  tend  to  strongly  deviate  from  the  properties  of  the  ideal 
vacuum.  The  transmissivity  of  the  radome  tends  to  decrease  as  the 
dielectric  constants  of  the  construction  materials  increase.  In  order 
to  meet  strength  requirements,  while  minimizing  the  perturbation 
caused  to  the  electromagnetic  field,  layer  construction  of  radomes 
has    been    developed. 

Lamination  construction  methods  allow  for  an  alternating 
combination  of  thin  layers  of  high  density,  high  strength  materials  of 
high  dielectric  constants,  with  thicker  layers  of  low  density,  lower 
strength  materials  of  low  dielectric  constants.  Electromagnetic  design 
of  laminated  radomes  is  more  complex  than  the  design  of  monolithic 
radomes  in  that  there  is  a  greater  number  of  boundary  conditions  to 
be  considered  in  establishing  the  lamination  layer  thicknesses.  The 
EMRAD  program  developed  by  Morgan  [Ref.  3]  currently  handles 
layered  dielectrics  with   up  to  five  layers. 

Figure  2.3,  adapted  from  the  Radome  Engineering  Handbook 
[Ref.  l:p.  46],  shows  several  types  of  laminated  wall  construction  that 
exhibit  the  alternation  of  layers  of  low  dielectric  constant  and  high 
dielectric  constant  materials.  The  two-ply  sandwich  is  the  simplest 
form  of  lamination  construction  with  a  thin  skin  of  strong,  dense 
material   of  high   dielectric   constant,   supported   internally   by   a   thicker 
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layer  of  porous  material  of  low  dielectric  constant.  This  design  is 
especially  attractive  at  low  frequencies  as  the  lateral  dimension  of 
the  thin  layer  of  dense  material  is  only  a  fraction  of  the  longer 
wavelength. 

The  A-sandwich  and  B-sandwich  are  two  variations  of  a 
stronger  three-layer  construction.  Although  the  A-sandwich  is 
superior  electrically  to  the  B-sandwich,  while  the  placement  of  the 
high  density  material  in  the  center  of  the  layer  construction  will 
strengthen  the  radome  wall  by  making  it  stiffer,  the  porous,  low 
density  material  on  the  exterior  of  the  radome  wall  will  not  sustain 
the   aerodynamic    forces   during   flight. 

The  B-sandwich's  inner  and  outer  skin  of  high  density  material 
and  the  core  of  low  dielectric  material  provide  stiffness  and  strength 
to  the  radome  structure  while  providing  a  dense  outer  skin.  This 
construction  has  a  higher  strength-to-weight  ratio  than  a  monolithic 
construction.  The  core  thickness  is  determined  by  the  strength 
requirements  and  provides  spacing  between  the  inner  and  outer 
skin.  Optimal  electrical  spacing  between  the  skins  would  by 
frequency  dependent  at  a  core  thickness  of  a  quarter  wavelength 
[Ref  l:p.  45]. 

The  C-sandwich  and  multilayer  radome  wall  constructions 
provide  for  very   high   strength   levels   which   may   be  required   at  high 
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Figure   2.3      Laminated   radome   wall   construction 
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operating  speeds  or  for  use  in  inclement  weather.  The  layer  material 
may  consist  of  plastic  laminates  and  ceramics.  Multilayered 
construction  is  attractive  where  extreme  structural  rigidity  or 
broadband   capabilities   are  primary   requirements.    [Ref.    l:p.   46] 

The  radome  interacts  with  the  electromagnetic  field  as  a  lens, 
modifying  the  field  as  it  passes  through  the  radome  structure.  At 
higher  frequencies,  this  interaction  can  predicted  by  the  use  of 
ray-tracing  techniques,  similar  to  optical  lens  analysis  [Refs.  2  and 
8].  Studies  have  been  conducted  at  higher  frequencies  to  design  a 
layer  of  dielectric  material  within  the  radome's  interior  structure  to 
act  as  a  correcting  lens  [Ref.  8]  as  shown  in  Fig.  2.4  .  This  package  of 
programs  will  allow  similar  analysis  of  lens  construction  at  lower 
frequencies,  where  the  radome  and  lens  parameters  are  within  the 
resonance   region   of  the  electromagnetic   fields. 

D.  MATERIALS  FOR  RADOME  CONSTRUCTION 

The  choice  of  material  for  construction  of  a  radome  is  another 
engineering  design  decision.  As  material  science  engineers  continue 
to  develop  new  materials,  particularly  low-weight,  high-strength 
ceramics,  the  range  of  materials  available  to  the  radome  designer 
increases.  The  general  properties  of  the  three  groups  of  materials 
which  are  in  common  use  are  listed  in  Ref.  9,  and  is  used  here  as  a 
source  for  material  selection.  The  groups  of  materials  discussed  are 
fiber-reinforced  organic  resins,  refractory  oxides,  and  glass  ceramics. 
Physical  properties  are  listed  in  Table  2.1. 
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Figure   2.4     Interior  radome  lens  of  dielectric   material 
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TABLE  2.1 


RADOME  MATERIAL  PROPERTIES 


FIBER 

REINFORCED 

GLASS 

REFRACTORY 

RESINS 

CERAMICS 

OXIDES 

Relative 

permittivity 

4-6 

5.5 

8-9 

(O 

Loss  tangent 

0.01 

0.003 

0.002 

(tan6) 

Density 


N/A 


160  lb/ft3  220  lb/ft3 


Ultimate 
Tensile 
Strength 
(UTS) 


40kpsi 


7  kpsi 


27  kpsi 


Compressive 
Strength 


20-40  kpsi  35  kpsi 


280  kpsi 


Maximum 
Temperature 


260  °C 


500  °C 


1700  °C 
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In  fiber-reinforced  organic  resins,  a  variety  of  resin  types  are 
used  to  impregnate  fiberglass.  Organic  resins,  such  as  phenolic, 
polyester,  and  epoxy  are  used  with  the  fiber  to  construct  the  radome. 
They  are  then  thermally  treated  to  set  the  structure.  Several 
manufacturing  techniques  are  discussed  in  Ref.  1.  Control  of  layer 
thickness  is  difficult  to  maintain  and  requires  finishing  to  avoid 
aberration.  Radomes  constructed  of  these  materials  have  poor 
resistance  to  rain  erosion,  but  are  relatively  cheap  and  easy  to 
manufacture,  even  in  large  sizes.  Maximum  temperatures  that  the 
material  can  sustain  limit  use  of  these  materials  to  missiles  with 
speed  of  less  than  Mach  1.5  to  Mach  2.0. 

Typical  glass  ceramics  used  in  radome  construction  are 
Pyrocream  and  glass-bonded  mica.  Compared  to  organic  radome 
materials,  ceramic  materials  have  much  higher  resistance  to  rain 
erosion  at  supersonic  speeds.  Ceramics  are  brittle  materials  and  the 
controls  needed  during  manufacturing  to  ensure  structural  integrity 
of  the  radomes  makes  them  expensive  to  build   [Ref.    l:p.  242]. 

The  refractory  oxide  used  in  radome  construction  is  usually 
alumina  (A1203),  but  magnesia  (MgO)  and  fused  silica  (Si02)  are  also 
used.  Refractory  oxide  manufacturing  is  also  expensive  due  to  the 
high  temperatures  used  in  forming  these  materials  to  the  desired 
radome  shape.  These  materials  have  high  dielectric  constants,  that 
can  undergo  a  change  of  up  to  ten  percent  between  room 
temperature  and  1700  °C.  Some  of  the  methods  by  which  refractory 
oxide    radomes    are    manufactured    are    slip    casting,    isostatic    pressing, 
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and  flame  spraying  [Ref.  l:pp.  298-313].  Oxides  have  a  high 
resistance  to  rain  erosion  and  wall  thickness  can  be  controlled  to 
within  0.002  inches,  allowing  for  extreme  control  of  aberration. 
Although  the  dielectric  constant  of  these  materials  is  high,  they 
provide  the  highest  degree  of  thermal  protection  and  place  the  least 
restriction   on    missile    speed. 
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IH.  ELECTROMAGNETIC  FIELD  THEORY 

A.  BACKGROUND 

The  analysis  of  the  interaction  between  a  missile's  radome  and 
the  electromagnetic  signal  that  is  radiated  or  received  by  that 
missile's  radar  is  a  problem  similar  to  those  undergoing  intensive 
investigation  in  electromagnetic  interaction  with  other  physical 
medium.  Ongoing  study  of  medium  interactions  include  other 
military  weapons  and  weapons  platforms,  forested  areas,  bodies  of 
water,  rain,  and  biological  tissues  [Refs.  3,  5,  and  10].  The  study  of 
electromagnetic  wave  penetration  through  missile  radomes  allows 
for  a  special  approach  to  the  solution  of  the  problem  for  the  common 
case  of  rotational  symmetry  in  the  construction  of  the  radome.  In 
the  high-frequency  regime,  the  solution  may  be  found  as  described 
by  Hirsch  and  Grove  [Ref.  2],  with  application  of  an  asymptotic 
approach  using  the  geometrical  theory  of  diffraction  (GTD),  which 
incorporates  ray-tracing.  When  lower  frequencies  are  considered, 
and  the  interaction  problem  is  in  the  resonance  region,  a  complete 
solution  of  Maxwell's  equations  is  required.  The  asymptotic  GTD 
method    yields    unsatisfactory    results. 

The  solution  of  Maxwell's  equations  for  the  radome  used  within 
this  project  is  divided  into  three  regions  as  illustrated  in  Fig.  3.1.  The 
three  regions  are  the  homogenous  spherical  inner  radome  core 
(designated  Region  I);  the  transition  region,  where  a  semi-annular 
conformal   mesh   is   generated   for  field   analysis   (designated   Region   II); 
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Figure   3.1     Three  spatial  radome  regions  for 
electromagnetic  field  analysis 
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and  an  unbounded  exterior  scattering  region  (designated  Region  III). 
The  total  solutions  for  the  fields  in  the  core,  Region  I,  and  the 
scattered  fields  in  the  exterior  region,  Region  III,  are  produced  via 
separate  spherical  harmonic  expansions.  Using  the  unimoment 
method  (Appendix  A),  the  modes  in  these  expansions  are  applied  as 
boundary  conditions  on  the  inner  and  outer  spherical  surfaces 
bounding  the  mesh  and  the  resultant  fields  for  each  such  mode  are 
found  using  the  finite  element  method.  A  similar  analysis  is 
performed   for   the   incident  fields   on   the   outer  surface. 

Using  superposition,  the  finite  element  solutions  are  then 
combined  using  the  same  unknown  coefficients  as  in  the  original 
expressions.  The  combined  finite  element  solution  is  equated,  in  the 
least  squares  sense,  to  the  original  expressions  along  the  inner  and 
outer  surfaces.  Enforcement  of  this  equality  allows  solution  of  the 
expansion  coefficients  for  the  core  and  exterior  regions. 

The  exterior  expansion  coefficients  are  used  by  EMRAD  to 
assemble  the  scattered  fields  and  bistatic  radar  cross  sections  (see 
Refs.  3  and  4).  The  core  region  spherical  harmonic  expansion 
coefficients  are  used  to  construct  the  fields  that  penetrate  the 
radome  to  illuminate  the  enclosed  antenna.  The  details  of  how  these 
interior  region  fields  are  generated,  using  the  core  region  spherical 
harmonics,  is  the  topic  of  the  next  section.  An  overview  of  the 
theory  used  by  Morgan  in  EMRAD  is  included  in  Appendix  A. 
Further  information  is  available  in  Refs.  3,  4,  5,  6,  and  7. 
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B.  GENERATION  OF  THE  CORE  REGION  FIELDS 

In  this  section  we  will  consider  the  construction  of  the  core 
region  fields  using  the  expansion  coefficients  provided  by  EMRAD. 
These  fields  are  assembled  and  plotted  in  various  ways  over 
user-selected  planar  surfaces  by  CORFLD. 

Within  the  spherical  core,  the  source-free  phasor  Maxwell's 
curl   equations,   with   suppressed   eJ®1    time  variation,  are 

-VxE=jcouH  (3.1) 

and 

VxH=jcoeE   ,  (3.2) 

where  u.  and  e  are  the  respective  permeability  and  permittivity.  The 
\x  and  e  will  usually  be  those  of  free-space  for  the  case  of  the  core 
region   in  a  radome. 

In  a  homogenous,  source-free  spherical  region,  such  as  the 
radome  core,  it  can  be  shown  [Ref.  11]  that  the  vector 
electromagnetic  field  can  be  represented  by  the  combination  of  two 
types  of  spherical  harmonic  fields,  transverse  magnetic  (TM)  and 
transverse  electric  (TE).  Using  the  spherical  coordinates  of  Fig.  3.2, 
these  TM  and  TE  fields  have  the  respective  properties  that  Hr  =  0  and 

Er  =  0.  Furthermore,  the  TM  and  TE  vector  fields  can  be  generated 
from  radially  directed   vector  potentials,   such   that   for  the  TM   fields, 

A   =  Arr  (3.3) 

and 
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Figure     3.2       Spherical  coordinates  with  unit  vectors 
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F  =  Fr  r  , 

where   the   radial   components    satisfy 


(3.4) 


f  A        \ 
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(  A       N 
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Ar 


Vr    J 


=  0 


(3.5) 


and 


^  ~\ 


+k" 


Fr 


vr  J 


=  0  . 


(3.6) 


These    equations     are    in     the    form    of  the    complex     scalar 
Helmholtz    equation, 

V2¥  +k2¥  =  0   .  (3.7) 

In   spherical   coordinates,   the   Helmholtz  equation  becomes   [Ref.    11], 
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+k2¥=0  .       (3.8) 


To    find    a    solution    to    this    differential    equation,    the    method    of 
separation  of  variables  is  used.     Letting 

¥  =  R(r)  6(9)  0>((j))  ,  (3.9) 

the    separation   of  variables   procedure   yields   the   system   of  ordinary 
differential    equations: 
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(3.10) 
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(3.1  1) 


and 


do; 

d<|>2 


+m20  =  0 


(3.12) 


The  R  equation  is  related  to  Bessel's  equation,  and  its  solutions  are  in 
the  form  of  spherical  Bessel  functions.     Spherical  Bessel  functions,  bn  , 

are  related  to  ordinary  Bessel  functions,  Bn,  as  indicated. 


R(r)   =  bn(kr)   = 


V2kr   X+l 


K) 


(kr) 


(3.13) 


The  0  equation  is  related  to  Legendre's  equation  and  its  solutions  are 
the   associated   Legendre   functions, 

0(6)  =  L™(cos0)   ->    Pnm(cos0)  ,  Qj;(cos9)   ,  (3.14) 

where    Pnm(cos  6)   are   the   associated   Legendre   functions   of  the   first 

kind   and   Qnm(cos   6)    are   the    associated    Legendre    functions    of   the 

second  kind.  The  associated  Legendre  functions  of  the  second  kind 
tend  toward  infinity  at  points  along  the  z-axis  and  must  be  excluded 
for  solutions  which  include  this  region. 
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The  <J>    equation   has   harmonic   sinusoidal   type   solutions 
0>((t>)   =  h(m(|>)->   cos(m<l>),  sin(m<|)),  e±im<t> 
=  c1cos(m(t>)+  C2sin(m(|)) 

=  C3e+jm<t,+  C4e-jm<1,  .  (3.15) 

Together,    these    separated    solutions    form    the    product    solutions 
to  the  Helmholtz  equation  as 

¥m>n   =  bn(kr)   Pnm(cose)  h(m<|))   .  (3.16) 

Equations  of  this  type  form  the  elementary  spherical  harmonic  wave 
functions.  The  general  solution  to  the  Helmholtz  equation  is  formed 
by  summation  of  linear  combinations  of  these  elementary  functions, 
one  for  the  TM  fields  and  one  for  the  TE  fields.  Each  of  these 
solutions  will  be  of  the  form 

y(r,e,(j>)  =  XX  Mkr>  Pnm(cose)  h(m<>)  •  <3-17) 

m     n 

The  summation  over  all  possible     m     and     n     is  limited   by   the 
nature  of  the  Legendre  polynomial   of  the  first  kind   in  that  Pnm    =  0 

for  n  <  m.  The  summation  is  also  limited  by  the  restriction  that  n 
will  not  equal  zero  when  m  is  zero.  For  both  n  and  m  equal  to 
zero,  then  the  field  evaluated  is  zero,  therefore  this  solution  is  also 
excluded.  Then  m  is  zero  and  the  summation  over  n  will  start 
with     n     equal  to  one.     This  leads  to  the  scalar  field  expansion 
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oo  oo 


^(r,0,<|))=X  X(cmncosm<t>+dmnsinm<|>)  P™(cos9)   bn(kr).        (3.18) 

m=0  n=m 

As  the  Helmholtz  equation  under  consideration  is   given   as 


then  for  the  TM  case 


V2¥  +k2¥  =  0, 


¥  =  &-  , 


(3.19) 


(3.20) 


and  for  the  TE  case 


F 

vp    _      r 

r 


(3.2  1) 


The    spherical    harmonic    solutions    of    the    Helmholtz    equations    for 
these  sets  of  modes  are  then  of  the  form 


Ar 


=  rVmfl   =  [(kr)  bn(kr)]  Pnm(cos0)   h(m<J>) 


(3.22) 


In  this  renormalized  form  of  the  solution  to  the  scalar  Helmholtz 
equation,  the  Riccati  (or  Schelkunoff)  spherical  Bessel  functions  are 
present.  The  Riccati  spherical  Bessel  functions  relate  to  the  spherical 
Bessel  function  and  normal  Bessel  functions  as  follows: 


JnOcr)  =  [(kr)bn(kr)]  =  ^-    B^(kr)   . 


(3.23) 


Therefore,  in   the   spherical  core  region,  the   solutions   for  Ar  and  Fr  are 
of  the  form 
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Ar 


oo  oo 


=  £    I>mncosm<j>+dmnsinm<i>)  Pnm(cos0)   Jn(kr)   .     (3.24) 
m=0  n=m 
n*0 


When    expressed    in    terms    of    exponential    Fourier    series    in    <J>,   the 
solutions  are  of  the  form 


Ar 

[Fr  . 


oo  oo 


=  I     Iamn    Pnm(cos0)  Jn(kr)  eJ 


m<t> 


(3.25) 


m=  — oo  n=m 
n*0 


The    coefficients    within    the    two    summations    are    determined    by    the 
solution  of  the  scattering  problem  provided  by  EMRAD. 

The  spherical  coordinate  electric  and  magnetic  field 
components  can  now  be  found  for  the  TE  and  TM  cases.  For  the  TM 
case,    where 


Htm   =    Vx(Arr) 


and 


(3.26) 
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For  the  TE  case,  where 
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The  total  fields  are  the  sum  of  the  TE  and  TM  parts 
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(3.3  6) 
(3.37) 

(3.3  8) 
(3.3  9) 

(3.40) 
(3.4  1) 


31 


1 

jwe 


(  <i 

dT2      '  ^ 


^    = L_    ™L    +    J_    I    ^r 


rsin9    9<}>        jcoe   r    9r99 


1    9F, 


1 


E*  =    7  3^  +  ~ 


1 


a2A, 


r     96        jcoe  rsin6    9r9<j> 


'  92 

TTFr+k2pr 

9rz 


=  _j_  a\  +  j^  i   a2Fr 


rsinG     9<})        jco^i   r    9r98 
1    9A,  1  1         92Fr 


(3.42) 

(3.4  3) 

(3.44) 
(3.45) 

(3.46) 
(3.47) 


r     90         jcon   rsin0    9r9<{) 
The   solution   for  the  coefficients   in   the   Ar  and  Fr   expansions    is 

performed  by  the  finite  element  and  unimoment  solution  method 
explained  in  Appendix  A.  The  finite  element  method  is  used  by 
EMRAD  and  applied  to  Region  II  between  the  inner  core  and  outer 
exterior  region.  In  this  finite  element  region  there  is  the 
inhomogeneous    axially    symmetric    structure    of   the   radome. 

EMRAD  generates  solutions  for  incident  plane  waves  with 
Poynting  vectors  at  incident  angles  specified  by  the  program  user, 
such  that  the  Poynting  vector's  direction  is  specified  by  the  unit 
propagation  vector  p  =  sin  a  x  +  cos  a  z.  As  any  plane  wave  can  be 
decomposed  into  the  sum  of  two  orthogonal,  linearly  polarized  plane 
waves,  EMRAD  evaluates  two  plane  waves  for  each  incident  angle,  a 
TM    plane    wave    and    a  TE  plane  wave,  as  shown  in  Fig.  3.3.    The  TM 
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Figure  3.3       Incident  Poynting  vector  with  TE  and  TM  plane  waves 
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and  TE  plane  waves  are  orthogonal  to  the  Poynting  vector  and  to 
each  other.  The  TE  plane  wave  is  always  in  the  y  direction.  The  TM 
wave  is  in  the  x-z  plane  and  directed  as  determined  by  the 
expression  cos  ax  -  sin  a  z.  For  incident  plane  waves  specified  in 
this  manner,  the  solution  for  the  component  values  of  the  E  fields  in 
the  core  of  the  radome  can  be  determined  from  taylored  equations 
for  the  TM  and  the  TE  incident  plane  waves.  For  the  TM  incident 
wave,  the  core  E -field   component   values   are  determined   as: 


^Er  (k0r)z  m=0 


n=m 
n*0 


p„m(cose) 

bnn(n+l)JD(kr)  n\       J 
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(3.4  8) 
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and 
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(3.50) 


For   the   TE   incident   wave,    the   core   E -field    component    values    are 
determined    as: 

ET(rM)=  P"    ,2IxmSin^(D)S 
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and 
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With  the  an  and  bn   coefficients  supplied   by  EMRAD  for  each  pair  of 

plane  waves,  TM  and  TE,  at  each  user-specified  angle,  CORFLD 
generates  the  interior  phasor  field  components  using  the  equations 
above.  The  total  magnitude  of  the  complex  vector  field  at  any  point 
within  the  radome's  core  is  expressed  as  the  square  root  of  the  sum 
of  the  squares  of  the  complex  absolute  values  of  the  field 
components.  This  total  magnitude  is  computed  and  displayed  for 
each    incident    wave    by    CORFLD,    using    the    an   and  bn     coefficients 

supplied  by  EMRAD. 

In  order  to  determine  the  perturbation  of  the  electromagnetic 
field  due  to  the  presence  of  the  radome,  CORFLD  generates  the  values 
of  the  incident  plane  wave,  then  compares  this  to  the  core  region 
field  components  calculated  by  EMRAD  and  CORFLD.  The  formulas 
for  the  undisturbed  incident  TM  and  TE  planes  waves  are 

E™=E0   e"jkop,r   (xcosa-zsina)  (3.54) 

and 

ETC=E0  e-jkoP#r    y    .  (3.5  5) 

The  magnitude  factor,   E0  =  1.0,  is  set  by  EMRAD  for  both  TM  and  TE 

waves.  The  unit  propagation  vector,  p,  includes  the  direction  of 
propagation   for  the  incident  waves   and  is  designated   as 

p  =  (x  sina+z  cosa)  .  (3.5  6) 
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This  vector  is  dotted  with  the  position  vector,  r,  for  each  point  on  the 
plane  of  interest,  as  defined  as 

r  =  (xx    +  yy+    zz)  .  (3.57) 

With  the  dot  product  evaluated  in  the  exponential,  the  TM  and  TE 
waves   are  specified   as 

E™=E0  e-Jk°(xsina+zcosa)  (xcosa-zsina)  (3.5  8) 

and 

ETE=E0    e-Jko(*sina+zcosa)     y      _  (3  5  9) 

The  relative  magnitude  and  phase  of  the  computed  field  in  the 
direction  of  the  theoretical  field  are  calculated  and  displayed.  This 
"scaled  dot  product"     is  defined  as 

Ed0t  =  ^  ,  (3-6  0) 

where    E    is   the   computed    field,    ej    is   a   unit   vector   pointing   in   the 

same  direction  as  either  the  TM  or  TE  incident  field  specified  in 
Equations   (3.58)   and   (3.59),   and 

Ej    =E0    e-JM*sina+zcosa)  (3  61) 

The  remaining  field  term  that  is  computed  and  displayed  by  CORFLD 
is  the  error  term.  The  error  term  is  that  portion  of  the  computed 
field  which  is  perpendicular  to  the  incident  field  vector.  In 
normalized   form   this   becomes, 


^error 


E-(E.ej)ei 


Ei 
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(3.62) 


The  relative  magnitude  of  the  error  term  is  calculated  and  displayed 
by  CORFLD  as  a  measure  of  electromagnetic  disturbance  caused  by 
the  presence  of  the  radome.  This  built-in  theoretical  field 
comparison  is  used  in  Chapter  V  as  a  measure  of  modelling  accuracy 
by  comparing  the  generated  theoretical  fields  to  the  fields  computed 
and  assembled  by  EMRAD  and  CORFLD  for  an  ideal  radome 
constructed  of  air  (er  =  1). 


C.  CANTING  THE  PLANAR  ARRAY 

The  planar  antenna  within  a  radome  can  be  canted  in  any 
direction  allowed  by  the  gimbal  system  that  is  used  with  the 
antenna.  In  order  to  evaluate  the  field  components  across  a  canted 
planar  antenna  within  the  radome  core,  the  coordinates  of  each 
point  on  the  antenna  must  be  converted  from  the  local  planar 
coordinate  system  to  the  global  radome  core  coordinate  system. 
CORFLD  uses  a  transformation  from  local  planar  to  global  core 
coordinates,  assuming  a  two-axis  gimbal  system  where  the  antenna 
is  rotated  about  the  original  z-axis,  as  shown  in  Fig.  3.4,  followed  by  a 
rotation  about  the  new  y-axis,  as  shown  in  Fig.  3.5. 

The  transformation  matrix  is  formed  by  multiplying  two 
matrices,  one  for  each  axis  rotation  [Ref.  12].  The  transformation  for 
the  rotation  about  the  z-axis  is  provided  by  the  matrix 


(3.6  3) 
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Figure  3.4  Rotation  about  the  z-axis 
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Figure  3.5   Rotation  about  the  new  y-axis 
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The    transformation    for   the   rotation    about    the   y-axis    is    provided    by 
the    matrix 
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The  complete   transformation   is   given   by   the  product 


(3.64) 
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(3.6  5) 


This  transformation  matrix  is  used  by  CORFLD  to  determine  the 
global  radome  coordinates  from  the  local  antenna  coordinates  for  a 
gimbaled  antenna  rotated  by  p  degrees  about  the  z-axis,  then  y 
degrees  about  the  new  y-axis.  With  this  computational  tool  the  fields 
that  are  received  by  an  antenna  within  the  core  of  the  radome  can  be 
determine  for  any  orientation  of  the  antenna  and  for  any  computed 
incident    field. 
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IV.  PROGRAM  EXECUTION  AND  DATA  FLOW 

A.  THE  RADOME  MODEL  DESIGN  PROCESS 

The  design  process  developed  by  this  thesis  is  represented  by 
Fig.  4.1.  The  radome  design  procedure  is  an  iterative  process  of 
selecting  radome  size  and  material  parameters,  radome  shaping 
functions,  and  placement  of  the  axis  origin  in  order  to  allow  for  an 
optimal  mesh  formation  by  EMRAD. 

The  flow  of  data  starts  when  RADSHP  is  used  to  design  the 
radome.  RADSHP  generates  a  material  file  which  is  ready  for  use  by 
EMRAD  and  a  structure  data  file  which  is  in  the  format  required  by 
CURVE  DIGITIZER.  CURVE  DIGITIZER  is  used  to  visually  check  the 
radome  structure  as  designed  in  RADSHP.  After  verifying  the 
structure  data  file,  the  program  EMCADIN,  developed  by  Connolly 
[Ref.  4],  is  used  to  translate  the  data  file  into  the  form  required  by 
EMRAD.  After  data  translation,  the  program  EMESH,  developed  by 
Morgan  [Ref.  3],  is  used  to  check  that  the  mesh  to  be  generated  in 
EMRAD  is  optimally  placed.  The  placement  of  the  origin  in  RADSHP 
when  generating  the  structure  data  file  can  be  modified  by  a 
placement  factor.  With  the  material  and  structure  files  are  prepared, 
EMRAD  is  executed  to  generate  the  interior  coefficients  of  expansion 
for  use  by  CORFLD.  A  sample  annotated  execution  through  the  data 
flow  follows. 
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Figure  4.1     Radome  model  design  process 
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B.  RADSHP  EXECUTION 

RADSHP  is  a  menu-driven  program  that  accepts  data  from  the 
keyboard.  When  the  executable  version  of  RADSHP  is  available, 
RADSHP  is  invoked  by  typing 

>  RADSHP 

at  the  DOS  prompt.  The  program  will  load  itself,  then  provide 
introductory    information. 

WELCOME  TO  RADSHP 

THIS  PROGRAM  GENERATES  RADOME  MATERIAL  AND  SHAPE  DATA 
FILES  FOR  USE  WITH  EMRAD,  BASED  ON  USER  INPUT. 

RADSHP  prompts  for  data  file  names  for  the  material  and  structure 
output  files.  Sample  keyboard  inputs  are  provided.  User  inputs  are 
preceded  by    "==>". 

TWO  FILES  ARE  REQUIRED  TO  HOLD  PROGRAM  OUTPUT  DATA. 


PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED  TO 
HOLD  THE  PROGRAM  MATERIAL  PARAMETER  DATA.  THE  EXTENSION 
.MAT  WILL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.MAT 

==>  MATFIL 


PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED 
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TO  HOLD  THE  PROGRAM  OUTPUT  DATA.  THE  EXTENSION  .DAT  WILL 
WILL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.DAT 

==>  STRFIL 

RADSHP  will  allow  the  construction  of  radomes  with  up  to  five  layers 
of  construction.  At  present,  each  of  these  layers  must  completely 
enclose   the   previous   layered   region. 

ENTER  THE  NUMBER  OF  LAYERS  (.LE.  5) 

==>  2 

RADSHP  prompts  for  the  parameters  of  each  layer  of  construction. 

Each  layer  of  the  radome  will  have  material  constants, 
an  outer  height  of  the  radome  layer,  H, 
a  outer  base  radius,  R,  and 
a  shape  function,  &  shape  factor  parameters  (if  appl.), 

The  first  "layer"  will  be  the  inner  core  of  the  radome, 
normally  filled  with  air. 

Layer  No:     1 

This  is  the  inner  core  of  the  radome 


Please  enter  the  radome  layer  height,  H, 
and  the  outer  base  radius,  R,  for  this  layer. 

==>    1.25,0.50 
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RADSHP  provides  a  menu  of  the  shaping  functions  available  for 
radome  construction.  The  curved  parabola  is  not  normally  used  for 
outer  layers  of  the  radome,  but  is  included  as  an  option  for  inner 
layer    construction. 

The  following  shaping  functions  are  available, 

1.  General  Ogive 

2.  Tangent  Ogive 

3.  von  Karman 

4.  Power  Series 

5.  Parabola,  curved 

6.  Parabola,  capped  by  pointed  apex 

Indicate  desired  shaping  function  by  entering  a  number 
1,2, ...  or  6 


==>  2 


Each  shaping  function  requires  input  of  their  respective 
shaping  factors.  These  are  prompted  for  and  then  listed  for  each 
layer  as  the  radome  is  constructed.  The  tangent  ogive  and  the  von 
Karman  functions  have  shaping  factors  that  are  determined  by  the 
height  and  base  radius  of  the  radome  layer;  these  factors  are 
determined  and  listed.  The  factor  designation  follows  the  form 
adopted  in  the  equation  in   Chapter  II. 

The  tangent  ogive  shaping  function  determines  two 
factors  for  the  radome  construction.  These  are  calculated 
based  on  the  height  and  radius  of  the  radome  layer. 
These  factors  are  F4  =    6.50  and  F5  =    6.38 
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LAYER     SHAPE  HEIGHT       RADIUS    Fl  F2       F3  F4  F5  F6       F7 

1  TAN  Ogive      1.25  .50  .00        .00     .00       6.50      6.38    .00      .00 


Layer  No:     2 

This  is  the  outer  layer  of  the  radome 


Please  enter  the  radome  layer  height,  H, 
and  the  outer  base  radius,  R,  for  this  layer. 

==>    1.30,0.55 

The  following  shaping  functions  are  available, 

1 .  General  Ogive 

2.  Tangent  Ogive 

3.  von  Karman 

4.  Power  Series 

5.  Parabola,  curved 

6.  Parabola,  capped  by  pointed  apex 

Indicate  desired  shaping  function  by  entering  a  number 
1,2, ...  or  6 

==>  2 

The  tangent  ogive  shaping  function  determines  two 
factors  for  the  radome  construction.  These  are  calculated 
based  on  the  height  and  radius  of  the  radome  layer. 
These  factors  are  F4  =    6.42  and  F5  =    6.28 
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LAYER     SHAPE  HEIGHT       RADIUS    Fl  F2       F3  F4  F5  F6       F7 


1  TAN  Ogive      1.25  .50  .00        .00     .00       6.50     6.38    .00      .00 

2  TAN  Ogive      1.30  .55  .00       .00     .00       6.42     6.28      .00    .00 

After  all  the  structure  forming  data  is  entered,  RADSHP 
prompts  for  the  material  relative  permittivity  parameter  for  each 
layer  of  construction.  Layer  1  is  the  inner-most  region  and  the 
permittivity    will    normally    equal    unity. 

Layer  No:     1 

This  is  the  inner  core  of  the  radome 
Enter  Er  for  this  Layer: 

==>   1.0 

Layer  No:     2 

This  is  the  outer  layer  of  the  radome 
Enter  Er  for  this  Layer: 

==>  4.0 

With  the  structure  and  material  parameters  entered,  RADSHP 
prompts  for  headers  to  identify  the  radome  constructed.  This  header 
information  is  stored  in  the  material  data  file  and  will  follow  the  data 
files   throughout  the  execution   of  the  programs   in   this   package. 
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THE  HEADERS  ALLOW  THE  USER  TO  IDENTIFY  THIS 
SET  OF  DATA  FROM  ALL  OTHER  SETS. 

PLEASE  ENTER  HEADER  #1  (64  CHARACTERS  MAX) 

==>  This  is  a  test  data  set  for  RADSHP  example 

PLEASE  ENTER  HEADER  #2  (64  CHARACTERS  MAX) 
==>  Header  info  stays  with  data  set  throughout 


With  all  the  data  input  and  calculations  completed,  RADSHP 
reports  the  generation  of  the  material  and  structure  files  and 
terminates. 

Computations  and  output  are  completed,  the  data  file,  MATFTL.MAT 
holds  the  material  specifications  for  input  to  EMCAD. 

The  data  file,  STRFIL.DAT,  holds  the  radome  shape  data  for 
viewing  in  the  CAD  package,  CURVE  DIGITIZER. 

The  structure  file  required  for  input  to  EMCAD  can  be 
generated  by  the  executing  the  program  EMCADIN, 
using,  STRFIL.DAT,  as  input  to  EMCADIN, 

The  output  file  from  EMCADIN  will  be  in  proper  format 
for  input  to  EMCAD  for  the  structure  file. 

Stop  -  Program  terminated. 
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C.  RADOME  CONSTRUCTION  VERIFICATION. 

The  construction  of  the  radome  by  RADSHP  is  verified  by 
entering  the  CURVE  DIGITIZER  package  and  displaying  the  structure 
data  file  generated  by  RADSHP.  Figure  4.2  is  a  copy  of  the  CURVE 
DIGITIZER  display  of  the  radome  shape  created  in  the  example 
execution  above.  Upon  review  and  acceptance  of  the  structure  data 
file,  EMCADIN  is  used  to  translate  the  data  file  into  the  form  required 
by  EMRAD.    EMCADIN  is  invoked  by  typing 

>EMCADIN 
at  the   system  prompt.      The  program  will  load  itself  and  displays   the 
greeting   heading   and  prompts  for  the  structure  data  file   as  input. 

WELCOME  TO  EMCADIN 

THIS  PROGRAM  ACCEPTS  PROPERLY  FORMATTED  INPUT  DATA 
FROM  CURVE  DIGITIZER  AND  CONVERTS  IT  TO  A  FORM 
WHICH  CAN  BE  USED  BY  EMCAD. 

PLEASE  PRESS  ENTER  TO  CONTINUE. 

INPUT  FILENAME 

INPUTFN  IS  THE  INPUT  FILENAME  OF  THE  FILE  CONTAINING 
THE  DATA  OF  INTEREST  WHICH  WILL  BE  INTERPOLATED  ON  AND 
CONVERTED  TO  A  FORM  WHICH  CAN  SUBSEQUENTLY  BE  USED  BY 
EMCAD 

PLEASE  ENTER  THE  NAME  OF  THE  INPUT  FILE.  THE  EXTENSION 
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Figure  4.2     Tangent  ogive  created  by  RADFLD 
as  displayed  by  CURVE  DIGITIZER 
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OF  THE  FILENAME  MUST  BE  .DAT,  I.E.  FILENAME.DAT 
==>  STRFIL.DAT 

EMCADIN  then  prompts  the   user  for  the  data  file  name  to  hold   the 
structure  file  that  is  prepared  for  EMRAD  execution. 

OUTPUT  FILENAME 

OUTPUTFN  IS  THE  OUTPUT  FILENAME  OF  THE  FILE  CONTAINING 
THE  DATA  WHICH  HAS  BEEN  INTERPOLATED  AND  CONVERTED  TO  A 
FORM  WHICH  WILL  BE  USED  BY  EMCAD. 

PLEASE  ENTER  THE  NAME  OF  THE  OUTPUT  FILE.  THE  EXTENSION 
OF  THE  FILENAME  MUST  BE  INCLUDED  I.E.  FILENAME.DAT 

==>  STRFIL.STR 

There  are  several  data  files  used  within  these  packages,  the 
default  file  name  extensions  should  be:  ".MAT"  for  material  data  files, 
".DAT"  for  structure  data  files  suitable  for  CURVE  DIGITIZER,  and 
".STR"  for  structure  files  suitable  for  EMRAD. 

After  receiving  the  data  file  names  to  work  with,  EMCADIN 
prompts  for  angular  resolution  requirements.  EMCADIN  uses  linear 
interpolation  to  provide  the  capability  to  generate  output  data  with 
better  resolution  than  the  input  data.  The  default  resolution  of  data 
generation  by  RADSHP  is  one  thousand  points  in  the  x-z  plane,  using 
180  points  of  angular  resolution  in  EMCADIN  is  sufficient  for  good 
results. 
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RESOLUTION 

DELTHE  IS  THE  USER  INPUT  VALUE  OF  THE  DESIRED  THETA 
RESOLUTION  IN  DEGREES. 

PLEASE  ENTER  THE  DESIRED  DELTA  THETA  VALUE  IN  DEGREES. 

==>  180. 

Upon  completion  of  the  data  file  translation  by  EMCADIN,  the 
data  files  required  for  execution  by  EMRAD  are  complete  and  in  the 
proper  format.  MATFIL.MAT  is  the  example  material  file  and 
STRFIL.STR  is  the  example  structure  file.  Before  executing  EMRAD, 
the  efficiency  of  mesh  generation  to  be  expected  during  EMRAD  can 
be  previewed  by  using  the  program  EMESH  developed  by  Morgan. 
Figure  4.3  shows  the  mesh  generated  by  EMESH  for  the  example 
radome,  as  displayed  by  CURVE  DIGITIZER. 

D.  EMRAD  EXECUTION 

With  the  data  files  generated  and  an  optimal  mesh  solution,  the 
data  flow  is  ready  to  initiate  EMRAD.  EMRAD  is  a  menu-driven 
program  developed  by  Connolly  and  Morgan  and  has  been  slightly 
modified  to  provide  the  interior  spherical  coefficients  of  expansion 
for  use  by  CORFLD.    EMRAD  can  be  invoked  by  typing 

>EMCAD 
at  the  DOS  prompt.     The  program  will  load  itself  and  provide 
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Figure  4.  3  Tangential  ogive  with  mesh  generated  by  EMESH 
as  displayed  by  CURVE  DIGITIZER 
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introductory    information. 

******  WELCOME  TO  EMRAD  ****** 

THIS  PROGRAM  COMPUTES  EM  FIELD  SCATTERING  FROM  PENETRABLE 
BODIES  OF  REVOLUTION  USING  THE  COUPLED  AZIMUTHAL  POTENTIAL 
(CAP)  FORMULATION  IN  CONJUNCTION  WITH  A  VARIATIONAL  FINTTE- 
ELEMENT  TECHNIQUE  AND  A  TO-REGIONAL  UNIMOMENT  METHOD.  THE 
NUMERICAL  SOLUTION  IS  PERFORMED  USING  A  TWO-SWEEP  RICCATI 
TRANSFORM  WITH  CIRCUMFERENTIAL  MARCHING. 
UPDATED  VERSION  OF  ORIGINAL  U.C.  BERKELEY  CODE. 
MODIFICATIONS  BY  M.A.  MORGAN         MAR  1987  -  DEC  1989 

E.M.  CONNOLLY      JUL  1987  -  APR  1988 
R.J.  VINCE  SEP  1989  -  DEC  1989  . 


PLEASE  PRESS  ANY  KEY  TO  CONTINUE. 


*******  EMRAD  INPUT  ******* 

EMRAD  ALLOWS  A  USER  TO  INPUT  THE  NECESSARY  INFORMATION 
ACCORDING  TO  THE  LEVEL  OF  HIS  EXPERTISE. 

1  NOVICE  LEVEL  -  GIVES  BRIEF  EXPLANATIONS/DESCRIPTIONS 

OF  INPUT  VALUES.  DESCRIBES  THE  REQUIRED 
FORMAT  FOR  THE  INPUT  VALUE.  GIVES 
TYPICAL  VALUES  WHERE  APPLICABLE. 

2  EXPERT  LEVEL  -  ASSUMES  THE  USER  IS  FAMILIAR  WITH  THE 

INPUT  PARAMETERS  AND  FORMATS.  SIMPLY 
PROMPTS  FOR  REQUIRED  INPUTS. 

3  EXIT  EMRAD 

54 


PLEASE  SELECT  THE  LEVEL  OF  EXPERTISE  BY  ENTERING  1,  2,  OR  3. 
==>   1 

The  names  of  the  data  files  generated  by  RADSHP  and  translated  by 
EMCADIN  are  entered  at  the  keyboard  when  prompted  by  EMRAD. 

EMRAD  REQUIRES  DATA  INPUT  FROM  TWO  FILES  PREPARED 

IN  ACCORDANCE  WITH  THE  USERS  MANUAL.  THESE  FILES  MUST 

CONTAIN  MATERIAL  PARAMETERS  AND  STRUCTURE  PARAMETERS, 

RESPECTIVELY. 

ENTER  MATERIAL  PARAMETER  FILE  (D:FILENAME.EXTENSION) : 

==>  MATFIL.MAT 

ENTER  STRUCTURE  DATA  FILE  (D:FILENAME.EXTENSION): 
==>  STRFIL.STR 

EMRAD  opens  both  files  and  displays  the  header  information  which 
RADSHP  placed  in  the  material  file. 

This  is  a  test  dataset  for  RADSHP  example 
Header  will  stay  with  dataset  throughout 

EMRAD  creates  several  datafiles.  The  first  will  be  the  diary  file  for 
the  execution  of  EMRAD.  This  file  will  have  the  ".OUT"  extension. 
The   second   data   file   generated   will   have   the   same  root   name   as   the 
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diary  file,  this  is  the  data  file  which  holds  the  spherical  coefficients  of 
expansion  for  the  core  field  generation.  This  second  file  will  have  the 
extension  ".INT"  and  will  be  used  as  input  to  CORFLD. 

THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT 
RESULTS  OF  ALL  NUMERICAL  CALCULATIONS  CONDUCTED  BY  EMRAD. 
THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO  EXTENSION  IS 
REQUIRED  OR  DESIRED.  EMRAD  AUTOMATICALLY  APPENDS  AN 
EXTENSION  OF  .OUT  TO  YOUR  FILENAME. 

A  SECOND  OUTPUT  DATA  FILE  WILL  BE  CREATED  TO  STORE  THE 
COEFFICIENTS  FOR  THE  INTERIOR  CORE  EXPANSION.  EMRAD  WILL 
APPEND  AN  EXTENSION  OF  .INT  FOR  THIS  DATA  FILE. 
PLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILES. 

==>  INTFIL 


EMRAD  was  developed  by  Morgan  and  Connolly  to  determine  the 
scattered  field  and  the  radar  cross  section  of  the  body  of  revolution 
under  study.  Four  output  data  files  are  created  for  input  to  separate 
graphing    routines. 

THE  OUTPUT  GRAPHICS  DATA  FILE  IS  THE  FILE  CONTAINING  THE 
OUTPUT  DATA  FROM  EMRAD  TO  BE  USED  AS  INPUT  TO  GRAPHING 
ROUTINES.  THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO 
EXTENSION  IS  REQUIRED  OR  DESIRED.   EMRAD  AUTOMATICALLY 
APPENDS  AN  EXTENSION  OF  .TMT,  .TMP,  .TET,  AND  .TEP  TO  YOUR 
FILENAME  AS  IT  PRODUCES  FOUR  OUTPUT  FILES  FOR  GRAPHICS. 

.TMT >  TM  INCIDENCE,  F-THETA 

.TMP >  TM  INCIDENCE,  F-PHI 
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.TET  — -->  TE  INCIDENCE,  F-THETA 
.TET  - — >  TE  INCIDENCE,  F-PHI 

PLEASE  ENTER  THE  FILE  NAME  OF  THE  OUTPUT  DATA  FILE. 

==>  SCATGRF 

EMRAD  prompts  for  a  graphics  caption  to  be  included  with  the  data 
set,  this  caption  is  passed  along  to  CORFLD  and  is  included  with  each 
plotted  figure.  This  caption  affords  the  user  the  ability  to  identify 
this  data  set  from  all  other  sets. 

THE  GRAPHICS  CAPTION  IS  A  PERSONALIZED  CAPTION  ALLOWING  THE 
USER  TO  IDENTIFY  THIS  SET  OF  GRAPHS  FROM  ALL  OTHER  SETS. 
THE  MAXIMUM  LENGTH  OF  THIS  CAPTION  IS  64  CHARACTERS. 
NOTE:  WHEN  USED  WITH  THE  GRAPHICS  PACKAGE,  THE  PROGRAM  IS 
ABLE  TO  DIFFERENTIATE  BETWEEN  UPPER  CASE  AND  LOWER  CASE 
CHARACTERS. 

PLEASE  ENTER  ANY  GRAPHICS  CAPTION  YOU  DESIRE 

==>  Graphics  caption 

As  discussed  in  Appendix  A,  EMRAD  utilizes  a  semi-annular 
conformal  mesh  in  order  to  conduct  electromagnetic  field 
calculations.  EMRAD  prompts  the  user  for  a  minimum  and  maximum 
mesh  density.  EMRAD  conducts  the  analysis  of  the  radome  structure 
using  a  linear  variation  of  density  from  maximum  density  decreasing 
to    minimum    density. 

DMIN  AND  DMAX  ARE  PARAMETERS  OF  MESH  DENSITY  IN  TERMS  OF 
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ELEMENTS/INTERIOR  LAMBDA.  INPUT  VALUES  ARE  EXPECTED  TO  BE 
REAL,  I.E.  THE  DECIMAL  POINT  MUST  BE  INCLUDED. 
TYPICAL  VALUES  ARE  DMIN  =  10.  AND  DMAX  =  15. 

PLEASE  INPUT  DMIN 

==>  15 

PLEASE  INPUT  DMAX 

==>  18 

EMRAD  allows  for  analysis  of  up  to  five  incident  fields.  EMRAD 
prompts  for  the  number  of  incident  angles,  and  then  for  the  value  of 
each  angle,  in  degrees.  EMRAD  uses  the  poynting  vector  direction  for 
the  incident  angle.  This  translates  to  an  angle  of  180  degrees  for  an 
angle  incident  at  the  radome's  appex,  135  degrees  for  an  angle 
incident  45  degrees  off  the  apex,  and  90  degrees  for  an  incident 
angle  off  the  beam  of  the  radome. 

THE  NUMBER  OF  INCIDENT  FIELD  ANGLES  IS  THE  TOTAL  NUMBER  OF 
INCIDENT  FIELDS  THAT  IMPINGE  ON  THE  OBJECT  OF  INTEREST. 
THIS  PROGRAM  ALLOWS  A  MAXIMUM  OF  FIVE  INCIDENT  FIELD 
ANGLES.  WHEN  ENTERING  YOUR  ANSWER  PLEASE  DO  NOT  INCLUDE  A 
DECIMAL  BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT,  I.E. 

NA  =  3  ,  NA  =  1 

PLEASE  INPUT  THE  NUMBER  OF  INCIDENT  FIELD  ANGLES. 

==>  2 
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ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  1 
==>  180 

ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  2 
==>  135 

EMRAD  performs  scattering  calculations  by  starting  at  6=0,  and 
progressing  clockwise  to  0  =  180  in  equal  increments,  the  menu 
allows  the  user  to  determine  the  size  of  the  increments. 

THE  NUMBER  OF  SCATTERING  FIELD  THETA  POINTS  DETERMINES  THE 
SPACING  BETWEEN  THETA  POINTS  DURING  EMRAD  ITERATIONS  AND 
CALCULATIONS. 

DELTA  THETA  =  180/  (NUMBER  THETA  POINTS  -  1)  SO... 

NUMBER  THETA  POINTS  =  37  >  DELTA  THETA  =  5  DEGREES 

NUMBER  THETA  POINTS  =  19  >  DELTA  THETA  =  10  DEGREES 

WHEN  ENTERING  YOUR  ANSWER,  PLEASE  DO  NOT  INCLUDE  A  DECIMAL 
BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT.  I.E.  NT  =  19 

PLEASE  INPUT  THE  NUMBER  OF  SCATTERING  FIELD  THETA  POINTS 
==>  19 

As  part  of  the  scattering  analysis,  EMRAD  allows  for  eight  angles  of 
user    observation. 
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THE  NUMBER  OF  PHI  ANGLES  IS  THE  TOTAL  NUMBER  OF  PHI  ANGLES. 
THIS  PROGRAM  ALLOWS  A  MAXIMUM  OF  THREE  PHI  ANGLES.  WHEN 
ENTERING  YOUR  ANSWER  PLEASE  DO  NOT  INCLUDE  A  DECIMAL 
BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT.  I.E. 

NP  =  3  ,  NP  =  1 

PLEASE  INPUT  THE  NUMBER  OF  PHI  ANGLES. 

==>  0 

EMRAD  truncates  the  modal  analysis  of  the  electromagnetic  field 
infinite  fourier  expansions  for  the  incident  plane  waves  and  the 
interior  core  as  well  as  the  exterior  scattered  fields.  Mathematical 
estimates  of  truncations  are  provided  and  the  user  is  prompted  for 
truncation  limits.  Higher  truncations  may  increase  the  modelling 
accuracy. 

ENTER  MSTOP  (.LE.  13) 
ESTIMATED  "MINIMUM"  VALUE  IS:       4 
==>  5 

Modal  truncations  for  the  number  of  spherical  harmonic  modes  in 
the  spherical  radome  core  and  external  scattering  region  are 
estimated.  The  user  is  prompted  for  truncation  limits  for  the 
spherical    harmonics. 
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ENTERING  NO.  OF  INTERNAL  AND  EXTERNAL  EXPANSION  MODES 
ESTIMATED  (KI*RMIN)  "MINIMUM"  INTERNAL  Nl  =     4 
ENTER  Nl  (.GE.  MSTOP): 
==>  5 

ESTIMATED  (KO*RMIN)  "MINIMUM"  EXTERNAL  N2  =     9 

ENTER  N2  (.GE.  MSTOP): 
==>  9 

Check  that  N1+N2  is  .LE.  35  ....  Otherwise  Abort 
==>  [RETURN] 

When  executing  the  numerical  solution,  EMRAD  requires  a  large 
amount  of  temporary  disk  storage  space.  The  amount  of  storage 
required  is  estimated  and  provided  to  the  user.  The  user  must  select 
a  disk  drive  that  has  adequate  memory  available  for  usage.  The  use 
of  a  ram  drive  configuration  can  increase  the  speed  of  execution  by 
EMRAD. 

*****  Temporary  Storage  of  Riccati  Matrices  ***** 
Estimated  Disk  Space  Needed  in  MegaBytes:    1 .930000E-KX) 
ENTER  Hard  Drive  or  RAM  Disk  Letter  (  C,  D,  E  ... ): 

==>  E 
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With  all  of  the  data  and  parameters  specified,  EMRAD  begins  the 
numerical  analysis  of  the  electromagnetic  interaction  of  the  radome 
structure  and  the  specified  incident  plane  waves.  Execution  is 
monitored   by   reports   to   the   terminal   screen. 

INDEXING  PROGRAM  THROUGH  VALUES  OF  M 

M-LOOP  ....  M  =      0 

CALL  MESH 

ENTER  I-LOOP;    NO  STEPS  TO  COMPLETE:      38 

1=       1        NO  STEPS  TO  GO:      37 

CALLLODER 

CALL  MARCH 

1=      2       NO  STEPS  TO  GO:      36 

EX  M-LOOP,  DATOUT 

***************  EMRAD  COMPLETED  ************** 

Stop  -  Program  terminated. 

Upon  termination  of  the  program  EMRAD,  the  coefficients  for  analysis 
of  the  radome  core  fields  have  been  stored  in  the  data  file,  ready  for 
CORFLD  execution. 

E.  CORFLD  EXECUTION 

CORFLD   is   a  menu  driven  program  that  accepts   data  from   the 
data    file    created    by    EMRAD    and    from    the    keyboard.       When    the 
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executable  version  of  CORFLD  is  available,  CORFLD  is  invoked  by 
typing 

>  CORFLD 

at  the  DOS  prompt.  The  program  will  load  itself,  then  provide 
introductory    information. 

WELCOME  TO  CORFLD 

CORFLD  ASSEMBLES  THE  FIELD  ON  THE  INTERIOR  OF  A  RADOME 
WITH  THE  COEFFICIENTS  SUPPLIED  BY  AN  EMRAD  OUTPUT  FILE. 

CORFLD  prompts  for  the  name  of  the  data  file  generated  by  EMRAD. 
Only  the  root  file  name  is  required,  CORFLD  attaches  the  ".INT"  to  the 
root  file  name. 

THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT 

RESULTS  OF  INTERIOR  FIELD  CALCULATIONS  CONDUCTED  BY  EMRAD. 

THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY. 

NO  EXTENSION  IS  REQUIRED  OR  DESIRED. 

CORFLD  AUTOMATICALLY  APPENDS  AN  EXTENSION  OF  .INT  TO  YOUR 

FILENAME. 

PLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILE. 

==>INTFIL 

Header  information  that  has  followed  the  data  flow  is  obtained  from 
the  ".INT"  data  file  by  CORFLD  and  displayed  on  the  screen.     Notice 
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that  the  graphics  caption  entered  during  EMRAD  execution  as  well  as 
the   headers   entered   during   RADSHP   execution   are   displayed. 

TEST  DATA  FOR  USER  EXAMPLE 

THIS  HEADING  FOLLOWS  THE  DATA  THROUGHOUT 

Graphics  caption 

CORFLD  reads  the  material  parameters  stored  in  the  ".INT"  file  and 
displays  them  on  the   screen. 

SRE=  ( 1 .000000,-5.000000E-07) 

SRU=  (1.000000,-5.000000E-07) 

KR  =  (1.000000,-1.000000E-06) 

CORFLD  reads  the  number  of  incident  angles  evaluated  by  EMRAD 
and  displays  each  angle  as  the  poynting  vector  angle  of  incidence  and 
also  as  the  angle  of  incidence  to  the  radome  exterior.  The  radome 
incident  angle  is  displayed  so  that  the  user  can  use  the  angle  as  a 
reference   when    canting    the   planar   antenna. 

Incoming  Inc.  Angle  #  1 

EMRAD  Poynting  angle  of  incidence  =       1 80.000000 
Radome  incident  angle  =    0.000000E+00 

Incoming  Inc.  Angle  #  2 

EMRAD  Poynting  angle  of  incidence  =      1 35.000000 
Radome  incident  angle  =       45.000000 

After  displaying  the  angles  of  incidence  available  for  consideration, 
EMRAD   prompts   the   user  for  the   angle  of  interest.      Only   one   angle 
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may  be  selected  for  each  execution  of  CORFLD.     In  order  to  consider 
other  angles,  CORFLD  must  be  re-executed  by  the  user. 

Select  Angle  of  Interest  By  Number 
==>  1 
The  angle  of  consideration  is  angle  #  1  =    0.000000E+00  degrees 

After  an  angle  is  selected  for  analysis,  the  user  may  enable  a  rotation 
of  the  planar  array. 

The  planar  array  is  set  by  default  to  lie  in  the 
x-y  plane,  with  boresight  in  the  z-direction. 

Do  you  want  to  change  the  array  boresight  ? 
(  "Y"  or  "N"  ) 

==>  N 

CORFLD  reads  the  data  from  the  ".INT"  file  and  assembles  the  interior 
field    components. 

MSTART=  0  MSTOP=  5 

Reading  Coefficients  and  Constructing 
Fields  for  Each  Azimuthal  "m"  Mode 

Terminating  M-Loop  at  m=l  for  +/-  Z-Axis  Incidence 

Azimuthal  Mode  "m"  =  0 
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Reading  Coefficients 

Azimuthal  Mode  "m"  =  1 

Reading  Coefficients 
Generate  Fields  Across  Array 
SPHERICAL  FIELDS  COMPUTED 

When  the  interior  fields  have  been  computed,  the  user  may  check  for 
the  field  component  values  for  any  point  on  the  planar  array. 

Check  Values  of  Coordinates  and  E-Field  ?  (Y/N): 

==>  Y 

Enter  Array  Triple  Index  to  Check  (J,L,I) 
J=  x-point  1,2,..,  50 

L=  y-point  1,2,..,  50 

1=  1  for  TM,  2  for  TE 

XP,YP,ZP:  -5.115592E-02    -5.115592E-02     0.000000E+00 
XC,YC,ZC:  -5.115592E-02    -5.115592E-02     0.000000E+00 

IE-RADI=    5.023 160E-01 
IE-THI  =    1.391270E-02 
IE-PHII=    5.022836E-01 
IE-TOTI=   7.104954E-01 

Look  at  Another  Point  ?  (Y/N): 

==>  N 

The   user  is   supplied   with   the   selection   menu   to   allow   access   to   any 
component    representation    of    the    interior    fields    across    the    planar 
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antenna  which  is  oriented  as  specified  by  the  user.  Component  phase 
and  magnitudes  of  the  interior  fields  as  well  as  the  total  field 
magnitudes    are    available, 

CORFLD  OUTPUT  GRAPHICS  PRESENTATION 
CORE  FIELD  REPRESENTATION  SELECTION  MENU 

Please  select  the  type  of  field  representation  of  interest 
CORE  will  display  the  following  field  representations: 

1.  TM  INCIDENT  FIELD,  E-THETA 

2.  TM  INCIDENT  FIELD,  E-PHI 

3.  TM  INCIDENT  FIELD,  E-RADIAL 

4.  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

5.  TE  INCIDENT  FIELD,  E-THETA 

6.  TE  INCIDENT  FIELD,  E-PHI 

7.  TE  INCIDENT  FIELD,  E-RADIAL 

8.  TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

9.  COMPARE  COMPUTED  FIELD  TO  INCIDENT  FIELD 

10.  CHANGE  ASPECT  RATIO 

0.    FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 
Indicate  your  selection  entering  "l","2",...,"10",or  "0" 
==>  4 

When  the  user  has  made  a  selection  from  the  menu,  the  three 
dimensional   (3-D)  plotting  routine  is  called  and   the  user  is   allowed   to 
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adjust    the    view    angle    for    best    observation    of    the    3-D    plot,    see 
Fig.  4.4, 

CALL  PLOTTING  ROUTINE 

A  black  and  white  'stick'  representation  as  well  as  a  color  3-D  fill 
representation  are  available  to  the  user.  The  'stick'  representation 
generates  the  clearest  output  for  printing  and  is  used  for  all  the  data 
presented    here. 

Enter  0  for  3D-stick  or  1  for  3D-fill  option 

==>  0 

(**  SEE  FIGURE  4.4  **) 

Default  view  angles  are:  phi=-45.,theta=70. 

ENTER  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE 

==>  9,9 

When  the  user  is  finished  with  an  electric  field  component,  the  new 
view  angle  "9,9"  is  entered  and  the  user  is  returned  to  the  core  field 
selection    menu. 

CORE  FIELD  REPRESENTATION  SELECTION  MENU 


9.    COMPARE  COMPUTED  FIELD  TO  INCIDENT  FIELD 
Indicate  your  selection  entering  "l","2",...,"10",or  "0" 
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The  datafile  used  to  generate  fields  is  intfil.INT 
TEST  DATA  FOR  USER  EXAMPLE 
THIS  HEADING  FOLLOWS  THE  DATA  THROUGHOUT 
Graphics  caption 

TM  INCIDENT  FIELD,     E-TOTAL  FIELD  MAGNITUDE 

1.58  f=- 


l.BB 


win.  Z  =  +.631;  wax.  Z  =  +.825;  average  Z  =  +.656 
Inc.  A  =  +.BB0;  Plane  G  =  +.888;  Uiew:  phi  =  +  45.8  theta=  +  75.8 


B8 


Figure  4.4    3-D  plot  of  E -total  field  magnitude  for  example  ogive  radome 
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When    the   user   wishes   to   see   the   amount   of  perturbation   caused    to 
the  incident  plane  wave  by  the  radome,  selection   "9"   is  entered, 


==>  9 


Selection  "9"  sends  the  user  into  Menu  2,  for  core  field  comparison  to 
the  ideal  theoretical  plane  wave  which  would  be  incident  on  the 
antenna    without    the   radome    present. 

CORE  FIELD  REPRESENTATION  SELECTION  MENU  2 

Please  select  the  type  of  field  representation  of  interest 
CORE  will  display  the  following  field  representations: 

1.  TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 

2.  TM  INCIDENT  FIELD,  Scaled  Error  component 

3.  TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 

4.  TE  INCIDENT  FIELD,  Scaled  Error  component 

Indicate  your  selection  entering  "1","2","3",  or  "4" 

==>  2 

CALL  PLOTTING  ROUTINE 

Enter  0  for  3D-stick  or  1  for  3D-fill  option 

==>  0 

Default  view  angles  are:  phi=-45.,theta=70. 

(**  SEE  FIGURE  4.5  **) 
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The  datafile  used  to  generate  fields  is  intfil.INT 
TEST  DATA  FOR  USER  EXAMPLE 
THIS  HEADING  FOLLOWS  THE  DATA  THROUGHOUT 
Graphics  caption 

TM  INCIDENT  FIELD,     Scaled  Error  component 

B.  5B  r- 


1.88 


win.  Z  =  +.813.  «ax.  Z  =  +.436;  average  Z  =  +.227 
Inc.  A  =  +.888;  Plane  G  =  +.888;  Uiew:  plu=  +  88. 8  theta=  +  78. 8 


Figure  4.5  3-D  plot  of  scaled  error  component  for  example  ogive  radome 
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ENTER  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE 
==>  9,9 


CORFLD  OUTPUT  GRAPHICS  PRESENTATION 
CORE  FIELD  REPRESENTATION  SELECTION  MENU 


0.    FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 

Indicate  your  selection  entering  "l","2",...,"10",or  "0" 
==>0 

*********  core  ANALYSIS  COMPLETED  ************** 
Stop  -  Program  terminated. 

When  the  user  is  finished  with  the  angle  of  incidence  selected,  "0"  is 
selected  at  the  selection  menu  and  the  program  execution  is 
terminated. 
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V.  PROGRAM  VALIDATION  AND  RESULTS 

In  order  to  validate  the  execution  of  EMRAD  and  CORFLD,  the 
tangent  ogive  radome  generated  in  Chapter  IV  was  modified  to  have 
a  relative  permittivity  equal  to  that  of  air,  er=1.0.  The  field 
coefficients  were  generated  by  EMRAD  and  the  core  fields  were 
computed  by  CORFLD  and  compared  to  the  unperturbed  incident  field 
as  listed  in  Chapter  III. 

Figures  5.1  and  5.2  show  the  magnitude  and  phase  of  the 
scaled  dot  product,  as  defined  in  Chapter  HI,  between  the  computed 
and  theoretical  values  for  the  TM  and  TE  incident  fields,  respectively, 
for  a  "nose-on"  incident  angle.  Each  comparison  shows  a  94.9% 
magnitude  matching  with  an  almost  constant  2.59  °  phase  differential 
across  the  array.  Figure  5.3  exhibits  the  scaled  error  component  for 
the  TM  and  TE  incident  fields,  with  an  average  magnitude  of  error 
component  of  0.4%.  The  comparisons  for  this  "nose-on"  angle  of 
incidence  show  that  there  is  good  agreement  between  the  theoretical 
and  calculated  field  components.  It  is  observed  that  the  TM  and  TE 
components  displayed  are  90  °  out  of  phase  with  each  other.  This  is 
expected  for  "nose-on"  incidence  since  the  TM  wave  has  its  direction 
component  in  the  -x  direction,  while  the  TE  wave's  direction 
component  is  in  the  y    direction. 

Figure  5.4  shows  the  magnitude  and  phase  of  the  scaled  dot 
product   between   the   computed   and    theoretical   values   for   the   TM 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

1.59  g- 

1.25  t- 


1.B8 


1.88 


win.  Z  =  +.98Z;  wax.  Z  =  +.996;  average  Z  =  +.949 
Inc.  A  =  +.BBB;  Plane  G  =  +.BBB;  Uiew:  phi=+45.B  theta=+75.  B 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

5.BB 


1.B8 


win.  Z  =  +1.73;  max.  Z  =  +2.95;  average  Z  =  +2.59 
Inc.  A  =  +.BBB;  Plane  G  =  +.BBB;  Uiew:  phi=+8B.  B  theta=+75.B 


Figure  5.1     Dot  product  magnitude  and  phase  for  TM 
computed  fields  determined  for  radome  composed  of  air 
compared  with  theoretical  plane  waves 

74 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

1.58  et 


1.80 


1.88 


win.  Z  =  +.982;  max.  Z  =  +.996;  average  Z  =  +.949 
Inc.  A  =  +.888;  Plane  G  =  +.888;  Uieu:  phi=+45.B  theta=  +  75.8 

The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

5.88  =- 


min.  Z  =  +1.73;  max.  Z  =  +2.95;  average  Z  =  +2.59 
Inc.  A  =  +.888;  Plane  G  =  +.888;  View:  phi  =  +  88. 8  theta=+75.8 


Figure  52     Dot  product  magnitude  and  phase  for  TE 
computed  fields  determined  for  radome  composed  of  air 
compared  with  theoretical  plane  waves 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


TE  INCIDENT  FIELD, 
8.18  i— 


Scaled  Error  component 


1.B8 


-8.85  ■— 


win.  Z  =  +.888;  wax.  Z  =  +.823;  average  Z  =  +.884 
Inc.  A  =  +.888;  Plane  G  =  +.888;  View!  phi=+18.  8  theta=+8B.  8 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


TM  INCIDENT  FIELD, 


Scaled  Error  component 
8.  IB  i— 


8.85    - 


1.88 


5« 


ttktotofe 


fiUw.     -«*it5 


m 


8.  5« 
1.88 


-8.85 


I 


win.  Z  =  +.888;  wax.  Z  =  +.823;  average  Z  =  +.884 
Inc.  A  =  +.888;  Plane  G  =  +.888;  Uiew:  phi  =+75.B  theta=  +  8B.  8 


Figure  5.3    Scaled  error  component  for  TM  and  TE 
computed  fields  determined  for  radome  composed  of  air 
compared  with  theoretical  plane  waves 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

3.88  |- 

2.58  t- 


1.88 


1.88 


win.  Z  =  +.548;  wax.  Z  =  +1.66;  average  Z  =  +.874 
Inc.  A  =  +68. 8;  Plane  G  =  +68.8;  Uie«:  phi=+45.B  theta=  +  75.B 

The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

78.88 


l.BB 


1.88 


min.  Z  =  -27.9;  max.  Z  =  +68.1;  average  Z  =  -2.84 
Inc.  A  =  +68.  8;  Plane  G  =  +68.8;  Uiew:  phi  =  +  45.B  theta=  +  75.B 

Figure  5.4    Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  air  compared  with  theoretical 
plane  waves  for  incident  angle  of  60  degrees 
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incident  field,  for  an  incident  angle  60°  from  the  "nose-on"  aspect. 
The  magnitude  of  the  scaled  dot  product  between  the  computed  and 
theoretical  TM  incident  wave  shows  modal  variation  around  the 
antenna  plane.  The  magnitude  of  the  computed  value  in  the 
direction  of  the  theoretical  plane  wave  varies  between  54.0%  and 
166%,  with  an  average  magnitude  of  87.4%.  The  phase  component  of 
the  dot  product  varies  across  the  plane  from  -27.9°  to  +60.1°,  with 
an  average  phase  of  -2.84°.  The  maximum  phase  differential  occurs 
at  the  edge  of  the  antenna  plane  closest  to  the  point  of  wave 
incidence,  as  the  incident  wave  is  defined  as  arriving  in  the  -x  -  z 
axis.  This  variation  seen  in  the  magnitude  and  phase  of  the  scaled 
dot  product  in  the  direction  of  the  theoretical,  undisturbed  plane 
wave  is  confirmed  in  the  scaled  error  component  of  Fig.  5.5  with  an 
amplitude  variation  from  8.4%  to  88.6%,  averaging  at  23.6%  of  the 
computed  field  evaluated  as  an  error  term  perpendicular  to  the 
theoretical    fields    expected. 

Figures  5.6  and  5.7  exhibit  the  behavior  of  the  calculated  TE 
60°  incident  plane  wave,  which  is  similar,  but  less  in  magnitude 
compared  to  the  calculated  TM  wave.  Less  modal  variation  is  seen 
along  the  exhibited  plane  of  the  antenna.  The  scaled  error 
component  has  a  symmetry  about  the  x-axis.  The  scaled  error 
amplitude  varies  form  0.3%  to  33.1%,  with  an  average  of  12.5%. 

The  form  of  the  error  components  as  analyzed  for  this  radome 
constructed  of  air  identifies  a  phase  error  in  the  determination  of  the 
computed  fields.  The  data  computations  of  the  EMRAD  subroutines 
have   been   previously   validated   by   Morgan.      These  validations 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


Til  INCIDENT  FIELD, 


Scaled  Error  component 
1.58 


l.BB 


1.  88 


min.    Z  =   +.084."    max.    Z  =    +.886;    average   Z  =   +.236 
Inc.    A  =   +68. 8;  Plane   G  =   +  68.  B;    Uieu:    plu  =  +  45.B  theta=  +  75.8 


Figure  5.5     Scaled  error  component  for  TM  computed  fields 
determined  for  radome  composed  of  air  compared  with 
theoretical  plane  waves  for  incident  angle  of  60  degrees 
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The  daiafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

2.  88  jr- 


1.88 


1.88 


win.  Z  =  +.868;  Max.  Z  =  +1.34;  average  Z  =  +.981 
Inc.  A  =  +68. 8;  Plane  G  =  +68.8;  View:  phi=+45.B  theta=  +  75.B 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

27.88  — 


-8.58 


-8.58 


1.88 


8.58 


8.88 


8.58 


1.88 


-33.88  l- 


nin.  Z  =  -18.2;  wax.  Z  =  +4.96;  average  Z  =  -.438 
Inc.  A  =  +68. 8;  Plane  G  =  +68.8;  View:  phi=+45.8  theta=  +  75.8 


Figure  5.6     Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  air  compared  with  theoretical 
plane  waves  for  incident  angle  of  60  degrees 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


TE  INCIDENT  FIELD, 
0.58 


Scaled  Error  component 


l.BB 


win.  Z  =  +  .BB3;  wax.  Z  =  +.331;  average  Z  =  +.125 
Inc.  A  =  +6B.B;  Plane  G  =  +6B.  B;  View:  phi=+ZB.  B  theta=  +  75.B 


Figure  5.7    Scaled  error  component  for  TM  computed  fields 
determined  for  radome  composed  of  air  compared  with 
theoretical  plane  waves  for  incident  angle  of  60  degrees 
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considered  the  magnitude  and  phase  of  the  scattered  fields  [Refs.  3 
and  4],  and  the  magnitude  of  the  core  fields  [Ref.  10]  along  the  main 
coordinate  axes.  This  work  is  the  first  to  plot  both  the  phase  and 
magnitude  variations  across  the  entire  two-dimensional  plane  and 
identifies  the  possible  need  for  increased  computational  accuracy  in 
the  numerical  routines  of  EMRAD.  There  is  also  the  possibility  of 
unresolved  programming  errors  in  CORFLD  which  should  be  further 
investigated. 

B.  COMPUTATION  RESULTS  FOR  VARIOUS  RADOME  MATERIALS 

EMRAD  and  CORFLD  executions  for  the  tangent  ogive  standard 
established  in  Chapter  IV  were  run  for  varying  material  permittivity 
according  to  the  material  properties  listed  in  Table  2.1.  The 
comparisons  of  the  computed  fields  to  the  unperturbed  theoretical 
incident  plane  waves  are  included  only  for  the  "nose-on"  incident 
angle  which  has  good  results  for  the  ideal  "air"  construction. 
Fiber-reinforced  glass,  with  er=4.0;  glass  ceramic,  with  er=5.5;  and 
refractory  oxide,  with  er=8.0,  are  included.  The  form  of  field 
comparisons  are  similar  in  shape  for  each  radome  considered,  only 
the  magnitude  of  the  plots  varies.  As  expected  from  the  discussion 
in  Chapter  II,  the  attenuation  of  the  incident  plane  wave,  received  at 
the  plane  of  the  antenna  in  the  core,  increases  with  increasing 
relative  permittivity.  When  the  radome  was  constructed  of  air,  the 
computed    magnitude    of    each    incident    wave    was    94.9%.       When 

fiber-reinforced   glass   construction   is  considered,   Figures   5.8,   5.9,   and 
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5.10,  the  average  magnitude  of  the  computed  wave  is  60.1%  and  the 
average  scaled  error  component  is  22.7%.  For  the  radome 
constructed  of  glass  ceramic,  Figures  5.11,  5.12,  and  5.13,  the  average 
magnitude  is  66.9%  and  the  average  scaled  error  component  is  13.3%. 
For  the  refractory  oxide  radome,  Figures  5.14,  5.15,  and  5.16,  the 
average  magnitude  of  received  field  is  87.6%  and  the  average  scaled 
error  component  is  20.6%. 
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The  datafile  used  to  generate  fields  Is  P4FIL. INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  fiber-reinforced  glass 

TM  INCIDENT  FIELD,  Conputed  Field  dot  Ideal  Field,  MAGN. 
1.5B  g- 
1.25  f- 
1.8B  |- 
B.  75  t- 


l.BB 


1.88 


win.  Z  =  +.475;  wax.  Z  =  +.811;  average  Z  =  +.681 
Inc.  A  =  +.BBB;  Plane  G  =  +.BBB;  View:  phi  =+35.B  theta=*7B. B 


The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  fiber-reinforced  glass 

TN  INCIDENT  FIELD,  Conputed  Field  dot  Ideal  Field,  PHASE 


8858804 


l.BB 


-56.88  - 


nrnin.  Z  =  -71.1;  wax.  Z  =  -27.4;  average  Z  =  -51.8 
Inc.  A  =  +.888;  Plane  G  =  +.BBB;  View:  phi=+18.B  theta=+85.B 


Figure  5.8     Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  fiber-reinforced  glass  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  p4fil.INI 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  fiber-reinforced  glass 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

1.58 


l.BB 


l.BB 


Min.  Z  =  +.475;  Max.  Z  =  +.811;  average  Z  =  +.6B1 
Inc.  A  =  +.BBB;  Plane  G  =  +.BBB;  View:  phi=  +  45.B  theta=+75.B 

The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 
Material  of  forwation  is  fiber-reinforced  glass 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


l.BB 


0-B.V 


l.BB 


win.  Z  =  -71.1;  Max.  Z  =  -27.4;  average  Z  =  -51.8 
Inc.  A  =  +.B8B;  Plane  G  =  ♦.BBB;  Vieu:  phi=+8B.B  theta=  +  8B.  B 

Figure  59     Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  fiber-reinforced  glass  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  fiber-reinforced  glass 


TM  INCIDENT  FIELD, 


Scaled  Error  component 


B.58  r- 


1.88 


win.  Z  =  +.813;  wax.  Z  =  +.436;  average  Z  =  +.227 
Inc.  A  =  +.888;  Plane  G  =  +.888;  View:  phi=+95.B  theta=  +  7B.  8 


The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  fiber-reinforced  glass 


TE  INCIDENT  FIELD, 
8.58  r- 


Scaled  Error  component 


l.BB 


win.  Z  =  +.813;  wax.  Z  =  +.436;  average  Z  =  +.227 
Inc.  A  =  +.BBB;  Plane  G  =  +.B88;  View:  phi=  +  18.  8  theta=  +  7B.  B 

Figure  510     Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  fiber-reinforced  glass  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =5.5 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  HAGN 

1.5B  e- 


1.88 


1.88 


min.  Z  =  +.668;  max.  Z  =  +.883;  average  Z  =  +.669 
Inc.  A  =  +.BB8;Plane  G  =  +.888;  View:  phi=+45.B  theta=+75.8 


The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =5.5 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

28.  BB 


l.BB 


B.5B   B.88  "^58 


-8. 58 


8 -00   B.58   r 


88 


-92. BB  "=- 
min.  Z  =  -77.5;  max.  Z  =  -47.5;  average  Z  =  -61.9 
Inc.  A  =  +.B88;Plane  G  =  +.BB8;  View:  phi=+45.B  theta=+75.B 


Figure  5.1  1     Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =  5.5 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN 

1.58  e- 


l.BB 


l.BB 


min.  Z  =  +.668;  max.  Z  =  +.883;  average  Z  =  +.669 
Inc.  A  =  +.B8B;Plane  G  =  +.888;  Uieu:  phi=+45.8  theta=+75.B 


The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =  5.5 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

28.88 


1.88 


8. 58   B.BB  -^58 


JL^J^B   8.58   r 


BB 


Inc.  A  =  +.BB8;Plane  G  =  +.BBB;  Vieu:  phi=+45.8  theta=+75.B 


Figure  5.12   Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =  5.5 

TM  INCIDENT  FIELD,     Scaled  Error  component 

8.58  r- 


1.88 


min.  Z  =  +.885;  max.  Z  =  +.312;  average  Z  =  +.133 
Inc.  A  =  +.BB8;Plane  G  =  +.888;  Uieu:  phi=+88.8  theta=+7B.B 


The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radone  is  tangential  ogive, 
constructed  of  ceramic,  Er  =  5.5 

TE  INCIDENT  FIELD,     Scaled  Error  component 
8.58  r- 


88 


min.  Z  =  +.885;  max.  Z  =  +.312;  average  Z  =  +.133 
Inc.  A  =  +.888;  Plane  G  =  +.888;  Uieu:  phi  =  +  18.8  theta=+7B.  8 

Figure  513     Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  refractory  oxide,  Er  =  8 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

2. 88  e- 


1.88 


1.88 


min.  Z  =  +.634;  Max.  Z  =  +1.B8.*  average  Z  =  +.876 
Inc.  A  =  +.BBB;Plane  G  =  +.888;  View:  phi=+45.B  theta=+75.B 


The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  refractory  oxide,  Er  =  8 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


1.88 


-83.88 

nin.  Z  =  -68.3;  max .  Z  =  -39.4;  average  Z  =  -52.4 
Inc.  A  =  +.BB8;Plane  G  =  +.B8B;  View:  phi=+88.B  theta=+75.8 


Figure  5.14     Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 
Mater ial  of  formation  is  refractory  oxide,  Er  =  8 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN . 

2. 88  f- 


l.BB 


8.5B 


1.88 


min.  Z  =  +.634;  max.  Z  =  +1.88;  average  Z  =  +.876 
Inc.  A  =  +.888; Plane  G  =  +.888;  Vieu:  phi =+55. 8  theta=+75.B 


The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  refractory  oxide,  Er  =  8 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


1.88 


8.58 


7.88 


^BlTea 


min.  Z  =  -68.3;  max.  Z  =  -39.4;  average  Z  =  -52.4 
Inc.  A  =  +.88B;Plane  G  =  +.888;  Vieu:  phi=+85.8  theta=+88.8 


Figure  5.15    Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  OXOFIL.IHT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  refractory  oxide,  Er  =  8 


TM  INCIDENT  FIELD. 


Scaled  Error  component 


1.50  m 
1.25 

i.ee 

B.75 
B.5B 


1.B0 


min.  Z  =  +.BB6;  wax.  Z  =  +.56B;  average  Z  =  +.2B6 
Inc.  A  =  +.BBB;Plane  G  =  +.BBB;  Uiew:  phi=+85.8  theta=+75.B 


The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  refractory  oxide.  Er  =  8 


TE  INCIDENT  FIELD, 

1.5B 

1.25 

1- 

l.BB 

B.75 

|- 

B.58 

Scaled  Error  component 


l.BB 


8.5B 


l.BB 


min.  Z  =  +.BB6;  max.  Z  =  +.56B;  average  Z  =  +.2B6 
Inc.  A  =  +.BBB;Plane  G  =  +.B8B;  Uieu:  phi=*25.B  theta=*75.B 

Figure  516     Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 


92 


VI.  CONCLUSIONS 

This  thesis  has  presented  a  useful  set  of  programs  to 
numerically  model  a  radome,  designed  of  specified  shape, 
construction,  and  material,  and  to  evaluate  the  perturbation  of  the 
electromagnetic  field  as  it  passes  through  the  modelled  radome  to 
the    antenna    within    the    radome    structure.  The    radome    design 

methodology  was  discussed  in  Chapter  II,  along  with  the  radome 
shaping  functions,  wall  construction,  and  choice  of  materials  for 
radome  construction.  The  program  RADSHP  is  introduced  and  the 
shaping  function,  as  well  as  the  shaping  factors  for  each  function  are 
explained.  Chapter  III  developed  the  electromagnetic  theory 
required  to  assemble  the  field  components  within  the  interior 
radome  core,  using  the  coefficients  provided  by  the  EMRAD  program. 
The  theory  required  for  transforming  the  coordinates  of  the  canted 
antenna  plane  to  the  core  coordinates  for  an  assumed  two-axis 
gimbal  was  also  discussed.  Chapter  IV  explained  the  data  flow  for 
the  radome  model  design  process  by  using  an  example  of  the 
program  executions  for  the  design  of  a  monolithic  tangent  ogive,  with 
computational  analysis  by  EMRAD  and  field  component  generation  by 
CORFLD.  The  verification  of  the  radome  design  and  optimal  placing  of 
the  core  origin  by  the  use  of  CURVE  DIGITIZER  software  is  also 
discussed.  Chapter  V  discussed  the  computational  errors  discovered 
in  the  3-D  plotting  of  the  field  components.     The   scaled  magnitude  of 
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the  computed  perturbed  plane  waves  for  the  "nose-on"  incident 
angle  agreed  within  95%  of  the  theoretical  value  and  the  phase  of  the 
computed  field  was  an  almost  constant  2.59°  out  of  phase  with  the 
theoretical  incident  plane  wave.  The  errors  increase  when  the  angle 
of  incidence  is  off  the  "nose-on"  aspect,  with  modal  amplitude 
variations  from  54%  to  166%  of  the  theoretical  plane  wave.  Further 
diagnostics  are  required  to  determine  the  source  of  this  error. 

This  is  only  the  first  stage  in  the  development  of  a 
mathematical  model  for  the  interaction  between  the  radome  and  the 
antenna  system  that  the  radome  is  designed  to  enclose.  The  model, 
in  its  present  form,  does  not  consider  the  presence  of  the  antenna 
structure  within  the  radome  core.  Further  work  is  needed  to 
determine  the  effect  of  the  antenna  structure  and  its  interaction  with 
the  radome.  When  such  a  model  is  developed  the  predicted  data  can 
be  compared  to  data  collected  by  a  physical  antenna  system.  The 
culmination  of  this  effort  could  be  the  analytic  development  of 
corrective  lenses,  as  referred  to  in  Chapter  II,  to  correct  or  minimize 
the  perturbation  of  the  incident  electromagnetic  field  caused  by  the 
radome    structure. 
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Appendix  A 

EMRAD  FIELD  THEORY 

The  program  EMRAD  uses  the  unimoment  method  mentioned  in 
Chapter  II,  as  developed  by  Morgan  and  Mei  [Ref.  5]  to  decouple  the 
interior  and  exterior  spherical  solutions.  The  expansion  coefficients 
for  the  electromagnetic  fields  within  the  radome  structure  are 
generated  as  part  of  the  scattering  solution.  CORFLD  uses  these 
coefficients  for  construction  of  the  interior  fields.  The  theory 
required  by  CORFLD  to  formulate  the  electromagnetic  fields  within 
the  radome  is  discussed  in  Chapter  II. 

The  code  developed  by  Morgan  for  EMRAD  uses  an  optimized 
variational  finite-element  algorithm,  in  conjunction  with  a 
tri-regional  unimoment  method,  to  provide  scattering  solutions  for 
each  of  multiple  incident  fields  impinging  on  the  specified  structure. 
EMRAD  represent  a  significant  improvement  in  computational  speed 
and  versatility  over  computer  codes  which  find  the  solution  to  the 
electromagnetic  scattering  and  penetration  problem  by  use  of  surface 
current  integral  equation  formulations.  This  advantage  results 
because  of  the  highly  sparse  matrices  embodied  in  the  finite  element 
solution  vis-a-vis  full  matrices  found  using  an  integral  equation 
formulation.  An  overview  of  the  finite  element  scattering  algorithm, 
along  with  the  concept  of  coupled  azimuthal  potentials  and  the 
unimoment  method  are  included  in  this  appendix.  Further 
information  is  available  in  references  3,  4,  5,  6,  and  7. 
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A.  COUPLED  AZIMUTHAL  POTENTIALS 

In  the  radome  construction  the  material  constitutive 
parameters,  e(r)  and  n(r),  are  invariant  to  the  azimuthal  coordinate  <J> 
due  to  the  axial  symmetry  of  the  radome.  This  axisymmetry  allows 
implementation  of  vector  field  generation  via  coupled  azimuthal 
potentials,  as  developed  by  Morgan,  Chang  and  Mei  [Ref.  7]. 

For  ease  of  solution  and  representation,  a  normalized 
cylindrical  coordinate  system  is  adopted,  as  defined  by  (R  ,  Z  ,  <|>)  = 
(k0  p,  kQ  z,  <|>)    where    (p,  z,  (J>)     are     standard     circular    cylindrical 

coordinates    and    k0  =  2n/X0      is    the   free-space    wavenumber   of   the 

time-harmonic  field.  The  equations  in  this  form  are  related  to  the 
standard  spherical  coordinates  (r,  6,0)  as  shown  in  Fig.  A.l  using  the 
substitutions 


R  =  k0    r   sine  (A.l) 


and 


Z  =  k0   r    cos6  .  (A. 2) 

The   total   electromagnetic   fields   are   decomposed  into   azimuthal 
modes  by  an  exponential  Fourier  series  in  the  $  -coordinate    such    that 

E(R,Z,<J>)  =  £   em(R,Z)  eJm*  (A.3) 

m=— oo 

and 

r,0H(R,Z,«>)  =  £    hm(R,Z)  e>m*  ,  (A.4) 


m=- 
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Figure  A.l    Rectangular    (x,  y,  z),  spherical  (r,8,  4>),  and 
normalized  cylindrical  (R,  Z,8)  coordinate  systems 
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where  r|0    =  120  n    Q.  is  used  to  normalize  the  magnetic  fields,  so  that 
h  =  T]0H. 

Using  the  axial  symmetry  of  the  radome  construction,  the 
modal  fields  decouple  when  the  field  expansions  are  substituted  into 
Maxwell's  equations,  yielding  [Ref.  7],  after  judicious  cross 
substitutions  the  following  modal  field  generating  equations  in  terms 
of  two  scalar  potentials,  TjC  R,Z,m  )  and  r2(  R,Z,m  ): 

<J)    xem   =  j   fn    (  m^xVIi    -  R^VFj    )  (A.5) 

♦  •  em      =     |.  (A.6) 
<t>xhm      =     j  fm  (   m<j)xVr2    +  RerVr]    )                     (A.7) 

♦  hm       =     £■  (A.8) 
where   the  two  dimensional   gradient   operator  is   defined   as 

V      =  [  R  (9/3r)  +  Z  (3/9Z)  ]  (A.9) 

and    fm   is  a  multiplicative  function  given  by 

im(R,Z)      =  [   ur(R,Z)  ^(R.29   R2-  m2    ]4      .  (A.l  0) 

Note  that  the  potentials,   1^  and  r2)   are  proportional   to   the   azimuthal 

components    of   the   electric    and   magnetic   fields,    em  and  hm.     These 

potential  functions  are  continuous  everywhere,  including  at  dielectric 
and  magnetic  interfaces.  This  property  of  uniform  continuity  in  the 
field  is  used  to  obtain  a  numerical  solution  of  the  fields. 
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The  total  fields  may  be  described  in  term  of  the  uniformly 
continuous  scalar  potentials  r^  R,Z,m  )  and  r2(  R,Z,m  )  such  that 

em(R,Z)  =  j  fm   (mVI]   -  Ru^xVr;,)   +  <J>  S-  (A.l  1) 

and 

hm(R,Z)  =  j  fm   (mVT2   +  Rer<>xVIi)  +  *  ^  .  (A.l  2) 

This  method  of  expanding  the  electromagnetic  fields  into  modal 
components  is  named  the  Coupled  Azimuthal  Potential  (CAP)  method 
due  to  the  fact  that  the  scalar  potentials  T^R,  Z,  m)  and  r2(R,  Z,  m) 
are  coupled  in  these  field  equations  and   are  not  separable. 

Due  to  the  symmetry  in  the  azimuthal  coordinate,  0,  in  radome 
design,  an  arbitrary  meridional  cross  section  of  the  radome 
construction  may  be  represented  in  an  (R,  Z)-plane  such  that  there 
exists  a  two-dimensional  solution  domain  S,  with  a  boundary  curve 
forming  the  surface  of  revolution  55.  It  has  been  shown  by  Morgan, 
Chang  and  Mei  [Ref.  7]  that  in  this  solution  domain  the  potentials,  Tj 
and  r2,  satisfy  the  following  formally  self-adjoint  pair  of  partial 
differential    equations    (PDE's), 

V-(fm    Re^lD-h-l^   +  mV-(f™4>xVr2)  =  0  (A.l  3) 

and 

R 

with  Dirichlet  boundary  conditions  (BC)  for  T}  and  T2  specified  on  bS 
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-mV-ft^xVIl)  +  V-(J™   Ru,  vr2)+i^r2   =  0   ,  (A.14) 


In  addition  to  this  PDE  formulation  there  is  an  easily  managed 
variational  principle  [Ref.  7]  established  using  either  an 
Euler-Lagrange  approach  or  a  stationary-theorem  approach.  The 
variational  principle  is  based  upon  the  existence  of  a  stationary 
functional  of  the  potentials,     Tj  and  r2,  and  their  gradients  such  that 

sf  =  jl  L(R,z;r1,r2,vr1,vr2)  dRdz  ,  (A.15) 

s 
where    the    Lagrangian    operator,    L,       is    unique    within    an    arbitrary 
constant   multiplier   and    arbitrary    independent   additive    function.      The 
Lagrangian   operator   has   the   quadratic   form 

L  =  ^   VTr     (RerVIi    +  m<|>xVr2) 
+      ^    VT2-    (RurVT2   -  mfxVIi) 

-     {  (er  Ii2   +   urT22    )  .  (A.16) 

With  this  analytical  approach  to  the  solution  of  the  electromagnetic 
field's  behavior  through  the  axially  symmetric  radome  structure,  the 
determination  of  the  field  on  a  point-to-point  basis  is  accomplished 
through    a    numerical    finite-element    algorithm. 

B.  FINITE  ELEMENT  ALGORITHM 

The    numerical    solution    of    the    electromagnetic    field    problem 

involving     scattering     by,     and     penetration     of,     the     radome     is 

accomplished     by    the    development    of    a    finite-element    algorithm 

using    the    mesh    geometry    in    Fig    A. 2.       This    configuration    in    the 

meridional     (R,  Z)-plane     uses     linear     triangular     elements     which  are 
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Figure  A. 2     Finite  Element  Mesh  Geometry 
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arranged  by  the  EMRAD  program  to  conform  to  the  surfaces  of  the 
axisymmetric  radome.  The  selection  of  surface  conformation 
simplifies  numerical  analysis  of  the  normal  derivative  discontinuities 
at  the  surface  boundaries.  The  two-dimensional  mesh  configuration 
that  is  developed  is  named  the  semiannular  conformal  (SAC)  mesh. 

Figure  A. 2  graphically  represents  the  separate  spatial  regions 
which  are  of  interest.  The  homogeneous  inner  spherical  core  of  the 
radome  occupies  Region  I.  The  inhomogeneous  axisymmetric  radome 
wall  construction  cross-section  is  in  Region  II.  The  unbounded 
exterior  region  exists  exterior  to  the  outer  sphere  in  the 
homogeneous  free-space  of  Region  III.  Note  that  the  surfaces  of  the 
inner  sphere,  Region  I,  and  the  exterior  Region  III,  bound  the  area  of 
interest  for  numerical  analysis  in  Region  II. 

The  electromagnetic  fields  within  the  homogeneous  inner  core 
of  the  radome  (Region  I)  and  the  unbounded  exterior  (Region  III)  are 
defined  in  terms  of  the  appropriate  spherical  harmonic  modal 
expansions  of  the  radial  TE  and  TM  field  components.  At  the 
surfaces,  r  =  a  and  r  =  b,  the  modal  elements  in  these  fields  are 
applied  as  sets  of  Dirichlet  boundary  conditions.  Then,  the  numerical 
solution  for  the  fields  in  the  finite-element  domain  of  Region  II  is 
generated  for  each  boundary  condition.  Finally,  the  unknown  modal 
expansion  coefficients  for  the  core  and  exterior  scattered  fields  are 
obtained  by  enforcing  continuity  of  the  electromagnetic  fields  at  the 
interfaces  along  r  =  t1  and  r  =  r2.     This  approach  to  numerical  solution 

of  the  fields  is  the  basis  of  the  unimoment  method.  In  the 
implementation    of    the    variational    principle    using    the    finite-element 
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method,  the  azimuthal  modal  field  components  are  selected  as  the 
fundamental     unknowns. 

In  Region  II  of  Fig.  A. 2,  the  SAC  mesh  is  composed  of  numerous 
interior  nodes.  Each  of  these  interior  nodes  is  surrounded  by  six 
triangular  regions  known  as  "elements".  The  azimuthal  field 
components  within  the  mesh  are  represented  as  basis  function 
expansions, 

N 

<D   •   em(R,Z)      =     X    em(n)  Un(R2)  (A.17) 

n=l 

and 

<j>   •  hm(R,Z)      =      X    hm(n)  un(R,Z)  .  (A.18) 

n=l 

Each  of  the  respective  sets  of  em(n)  and  hm(n)  represent  one  of  the  N 
complex  nodal  values  found  in  the  azimuthal  fields.  The  finite 
element    basis    functions    un(R,   Z)   are   area   coordinate   linear   pyramid 

functions  which  each  assume  a  value  of  unity  at  their  associated 
node,  n,  and  zero  at  each  of  the  neighboring  nodes.  There  is  an 
assumption  of  linear  variation  within  each  of  the  triangular  elements 
in   the   surrounding   group   such   that 

un(R,Z)  =  a^R  +  bMZ  +  cM   .  (A.19) 

The  real  constants  (a,  b,  c)  corresponding  to  each  node,  n,     in  element 

£     is  found  in  terms  of  the  nodal  (R,  Z)-coordinates  of  the  element. 

The     variational     finite-element     technique     is     initiated     by 

substituting    the   linear   shape-function   expansions   in   Equations   (A.17) 

and    (A.18)    into    the    stationary    Euler-Lagrange    functional,    SF,    as 
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described  in  Equation  (A.  15).  The  stationary  point  of  SF  is  then 
computed  by  imposing  the  null  derivative  condition  with  respect  to 
each  of  the  unknown  coordinates  representing  the  interior  nodal 
fields,  while  constraining  the  "known"  boundary  nodal  fields.  The 
resulting  linear  system  describing  the  nodal  fields  is 

I{em(i)JJ   Kn  RerVfRuJ.VtRUiJ-Re^uJdRdZ 

i=l  Sn 

+hm(i) jj  [m^VfRuJ-^xVtRuOldRdZHO  (A.20) 

and 

N 

I{hm(i)JJ   [^nRuTV(Run)V(Rui)-RuTunui]dRdZ 
i=i  sn 

-em(i)JJ  [m41V(Run)-0xV(Rui)]dRdZ}=Q  (A.2  1) 

Sn 

for  each  value  of  the  internal  nodes,  n,  where  i  spans  all  nodes, 
including  boundary  nodes.  As  these  above  equations  are  enforced 
for  each  internal  node,  they  involve  integrations  over  only  the 
elements  connected  to  this  node,  Sn,  and  therefore  relate  the  value  of 
em(n)  and  hm(n)     only  to  adjacent  nodal  values,  em(i)  and  hm(i).     Then, 

in  the  SAC  mesh  implementation,  in  the  case  of  each  internal  node 
and  the  six  adjacent  nodes,  the  above  equation  will  relate  14 
adjacent  node  azimuthal  field  values,  some  of  which  may  be  specified 
at    the    boundary. 


104 


The  computation  for  each  nodal  value,  as  indicated  by 
Equations  (A. 20)  and  (A. 21),  results  in  a  sparse  global  system  matrix 
which  contains  generally  complex  elements.  At  each  node  the  mean 
values  of  the  spatially  variable  er(R,  Z)  and  \iT(R,  Z)  are  used  as  the 
material  constants  in  that  triangular  element.  Applying  the  Dirichlet 
boundary  conditions  to  the  SAC  mesh  shown  in  Fig.  A. 2,  Morgan  and 
Mei  [Ref.  5]  found  the  numerical  finite-element  solution  for  the 
electromagnetic  fields  in  Region  II  using  a  two-sweep  block-by-block 
matrix  inversion  algorithm  related  to  the  Riccati  transform  for 
difference  equations.  At  each  constant  colatitude,  9  =  6k,  a  vector  of 
internal  nodal  azimuthal  field  values  along  the  ray  from  t1  to  r2  is 
defined  by  |3  k  .  The  resultant  global  system  matrix  defined  by 
Equations  (A. 20)  and  (A. 21)  has  a  diagonal  tri-block  submatrix 
structure,  which  can  be  represented  locally  by  relationships  between 
adjacent  radial   node  vectors  in   the  mesh, 

Ak    Pk_i   +  Bk    fc    +  Q    fc+1    =  ak   ,  (A.22) 

such  that  Ak,  Bk,  and  Ck  are  sparse  submatrices.  These  form  the 
tri-block  diagonal  structure  of  the  global  system  matrix  while  a  k  , 
the  local  driving  function,  is  generated  from  Equations  (A. 20)  and 
(A. 21)  using  the  known  Dirichlet  BC's  at  r=a  and  r=b. 

The  system  solution  in  the  two-dimensional  meridional  (R,  Z) 
plane  is  performed  in  a  step-wise  manner  beginning  at  the  positive 
Z-axis  (0k  =  0!  =0)  and  proceeding  in  equal  marching  increments  of 
A6.     At  each  0k  ,      the    SAC    mesh    is    constructed    using    a    minimum 

distortion    scheme    to    conform    to    the    approximate    piecewise    smooth 
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surfaces.  At  each  0k  increment,  as  the  mesh  is  constructed  in  this 
sequential  circumferential  manner,  the  three  block  matrices  are 
generated   which   relate   adjacent   radial   node   vectors. 

During  the  incremental  marching  of  9k  ,  the  block-by-block 
Riccati  matrix  inversion  is  implemented.  The  Riccati  matrix,  Rk  ,  and 
the  associated  vector,  Sk  ,  are  related  to  the  vector,  p  k  >  by  the 
transformation 

Pk-i   =  Rk  Pk  +  Sk    .  (A.2  3) 

Substituting  this  Riccati  relation  into  Equation  (A. 22),  solving  for  P  k 
and  returning  to  the  form  of  the  Riccati  transform  yields  the 
recurrence    relations 

Rk+1   =  -(  Bk   +  Ak.R,    r1.^  (A.24) 

and 

Sk+1   =  (  B,   +  Ak.Rk  )-!.(  pk    -  Ak.Sk   )  .  (A.25) 

Using  these  relations,  the  Riccati  transform  sweeps  forward 
from  initial  conditions,  using  the  equations  of  step-wise  evolution  in 
Equations   (A.24)  and   (A.25)   to  sequentially  generate  and   store  the   Rk 

matrices  and  the  Sk  vectors.  This  is  called  the  "first  sweep".  The 
solution  vectors  P  k  are  then  obtained  by  using  the  transform  in 
Equation  (A.23)  to  perform  a  "backsweep"  with  calculations  using  the 
stored  Riccati  arrays.  This  sequentially  localized  technique  for 
handling    the    Dirichlet    boundary    value    problem    blends    easily    with 
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the    unimoment    method    which    is    used    to    solve    the    electromagnetic 
field  in  Regions  I  and  III  of  Fig.  A. 2. 

C.  UNIMOMENT  METHOD 

The  unimoment  method  is  a  computational  procedure  which 
provides  for  the  numerical  solution  of  the  incident,  scattered  and 
spherical  core  electromagnetic  field  expansions  when  dealing  with 
the  constructed  SAC  finite-element  mesh  of  Fig.  A. 2  The  numerical 
solution   is   generated   by  representing   the   core   (r<r2)  and  the  exterior 

scattered     (r>r2)    fields    in    terms    of    radial    TE    and    TM    spherical 

harmonic    expansions,    each    with    unknown    core    and    scattered    modal 

field   coefficients.      In   order  to   solve   for   these   unknown   coefficients, 

the   core   expansion    of   Region   I    is   applied   as   a   boundary   condition 

along  r  =   a,  while  the  boundary  condition  at  r  =  b  is  determined   as 

the  addition   of  the  incident  field  and  the  scattered   field  expansion   in 

Region   III. 

Superposition     is    then     applied     to    the    boundary     conditions, 

composed   of  the  core  field  expansion   and   the  exterior  field  expansion 

in   order  to  find   a  complete   solution.      The  complete  field   solution   is 

found   from   the   addition   of  the   partial   field   solutions   generated   from 

two  sets  of  applied  BC's.     The  first  set  of  boundary  conditions  applied 

is  composed  of  the  modal  spherical  harmonic  core  fields  along  r  =  a 

with   zero  fields   along  r  =  b.     The  second   set  of  boundary   conditions 

contains    the    spherical    harmonic   exterior   scattered    and    incident   fields 

at  r  =  b  with  null  fields  at  r  =  a.     These  multiple  boundary  conditions 

generate     multiple     solutions     at    the     nodes     along     both       boundaries, 
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r  =  fj  and  r  =  t: 


The    unknown    core    and    scattered    modal    field 


coefficients  are  obtained  by  setting  the  weighted  superposition  of 
interior  solutions  equal  to  the  original  field  expansions.  This  equality 
is  enforced  along  the  boundaries  at       r  =  t{  and  r  =  r2  ,   utilizing   least 

squares. 

The  success  of  the  unimoment  method  is  contingent  upon  the 
capability  to  analytically  represent  the  electromagnetic  fields 
everywhere  outside  of  the  finite  element  mesh  region.  These 
representations,  as  developed  for  Regions  I  and  III,  are  given  as 
infinite  series  of  spherical  harmonics  which  must  be  truncated  to 
facilitate  numerical  calculations.  The  plane-wave  fields  that  are 
incident  on  the  radome  are  decomposed  in  spherical  coordinates  into 
azimuthal    modes    as 


M 


E'(tM)  =1    e'     (r.6) 


m=0 


jsin(m<j>) 
cos(m<|>) 


(A. 26) 


and 


M 


TloHVrAO)  =X    h'     (r,6) 


m=0 


COS(m(f)) 


_  jsin(m<J)) 


(A. 27) 


where       M       is    the       mathematical    truncation    index    of    the    series 

representation.       The    bracketed    sinusoids    represent    the    polarization 

convention  used  for  TE  and  TM  modes,  respectively. 

During    computations,    the    numerical    solution    of    the    unknown 

core    and    scattered    modal    field    coefficients    is    run    concurrently    with 

the    sequential    Riccati    transform   finite-element   solution   in   Region   II. 
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On  the  forward  sweep,  boundary  conditions  for  both  the  incident  and 
the  modal  fields  are  generated  only  at  the  local  node  row  being 
considered  in  that  particular  sweep  step.  On  the  backsweep  the 
nodal  values  for  each  applied  boundary  condition  are  obtained  along 
the  r  =  rj  and  r  =  r2    contours.      The   corresponding   incident   field    and 

spherical  harmonic  field  at  each  node  are  generated,  and  then  the 
nodal  residual  fields  are  formed.  During  the  backsweep,  inner 
products  integrations  are  accumulated  and  their  contributions  are 
added  to  a  field  moment  matrix  and  the  incident-field  driving  vector, 
which  corresponds  to  each  incident  field.  The  field  expansion 
coefficients  for  each  incident  field  are  obtained  by  multiplying  the 
inverse  of  the  field  moment  matrix  with  that  field's  driving  vector. 
With  the  field  expansion  coefficients  calculated  and  collected,  the 
electromagnetic  field  within  the  radome  core  can  be  determined  on  a 
point-to-point  basis   as  developed  in   Chapter  II   . 
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c 


Appendix  B 
RADSHP  PROGRAM  CODE 


************************************************************** 


PROGRAM  RADSHP 
C 
C 
C 

C  GENERATING  THE  MATERIAL  FILE  AND  DATA  STRUCTURE  FILE 
C        FOR  A  LAYERED  DIELECTRIC  RADOME 

C  MATERIAL  FILE  PREPARED  FOR  EMRAD 

C  DATAFILE  PREPARED  FOR  CURVE  DIGITIZER 
C        USE  EMCADIN  TO  GENERATE  STRUCTURE  FILE  FOR  EMRAD 

C  Written  by  R.J.Vince  Nov-Dec  89 
C 

Q  ************************************************************** 

C      VARIABLE  DECLARATIONS 

*************************************************************** 

REAL  X(5,1001),  Z(5,1001),  PI,  DELTHETA,  DEGTORAD 

REALTHETA,C(5),R(5),H(5),DELZ(5),SHFT^ACTORJ;(5,7) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESOLV,  RESPTS 

COMPLEX  ER(5),UR(5) 

CHARACTER*  10  SHAPE(6) 

CHARACTER*64  HDR1.HDR2 

CHARACTER* 8  YNAME 

CHARACTER*  12  MATFIL.DATFIL 

CHARACTER*  1YN 

SLARGE 

C 

Q  *************************************************************** 

C      INITIAL  VALUES 

Q  *************************************************************** 
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PI=3. 14 15927 
DEGTORAD  =  PI/180. 
SEPRATN  =  999990 
RESOLV=1000 
RESPTS=RESOLV+l 
NLMAX  =  5 

C     ITITIALEE  SHAPE  FUNCTION  FACTOR  ARRAY 
C 

D067I=1,NLMAX 
DO  65  J=  1,6 
F(U)=0.0 

65         CONTINUE 
67         CONTINUE 
C 

C 

Q  ***************+**********+****************»************+»******+* 

C      INPUT  NECESSARY  DATA 

c 

C     PROMPT  FOR  OUTPUT  FILENAMES 
C 

WRITE(MIO) 
1 10        FORMAT(//9X,'WELCOME  TO  RADSHP', 

1  ///9X,THIS  PROGRAM  GENERATES  RADOME  MATERIAL  AND  SHAPE  DATA', 

2  //7X,'FTLES  FOR  USE  WITH  EMRAD,  BASED  ON  USER  INPUT.', 

3  ////9X/TWO  FILES  ARE  REQUIRED  TO  HOLD  PROGRAM  OUPUT  DATA.', 

4  ///9X,'  PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED  TO  ', 

5  /7X/HOLD  THE  PROGRAM  MATERIAL  PARAMETER  DATA.  THE  EXTENSION'. 

6  /7X/.MAT  WILL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.MAT7) 
READ(VOO)  YNAME 

NCHAR  =  INDEX(YNAME,'  ■)  -1 
MATFTL  =  YNAME(1:NCHAR)  //  '.MAT 

OPEN(  1  JTLE=MATFTL,STATUS=XJNKNOWN',ACCESS='SEQUENTIAL', 
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1      FORM='FORMATTED) 
WRITE(Mll) 
1 1 1         F0RMAT(///9X,*PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED  ', 

1  /7X.TO  HOLD  THE  PROGRAM  OUTPUT  DATA.  THE  EXTENSION  .DAT  WILL' 

2  f7X,rWJLL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.DAT/7X) 
READ<M00)  YNAME 

NCHAR  =  indexcyname;  0  -1 

DATFIL  =  YNAME(1:NCHAR)  //  '.DAT 

OPEN(2,FlLE=DATFIL,STATUS='UNKNOWN'^CCESS='SEQUENTIAL') 
1      FORM=*FORMATTED') 


C 

C     PROMPT  FOR  #  OF  LAYERS  WITHIN  THE  RADOME 


C 


5  WRITE(*,114)NLMAX 

1 14  FORMAT(//7X,'ENTER  THE  NUMBER  OF  LAYERS  (LE.  ',I2;)7/7X) 
READ(*,*)  NL 

IF  (  (NL  XT.  1)  .OR.  (NL  .GT.  NLMAX) )  THEN 
WRITE(M  15)NLMAX+1 

1 15  FORMAT(//7X,'NUMBER  OF  LAYERS  IS  OUT  OF  RANGE,((k#<',I2/)',//7X) 
GOTO  5 

ENDIF 


C 

C     LOOP  TO  INPUT  LAYER  INFORMATION 


C 


WRITE(*,7) 

FORMAT(/8X,'Each  layer  of  the  radome  will  have  material  consts,', 

1  /lOX.'an  outer  height  of  the  radome  layer.H,', 

2  /10X,'a  outer  base  radius,R,  and, 

3  /lOX.'a  shape  function,  &  shape  factor  parameters  (if  appl.),', 

4  //lOX.'The  first  "layer"  will  be  the  inner  core  of  the  radome, ', 

5  /10X,'normally  filled  with  air.') 
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D012L=1,NL 
39         CONTINUE 

WRTTE(*,9)  L 
9  FORMAT(/,'l     Layer  No:  *J3) 

IF(L.EQ.1)THEN 

WRITE(*,*) '  This  is  the  inner  core  of  the  radome' 
ELSEJF  (LEQ.NL)  THEN 

WRITE(*,*) '  This  is  the  outer  layer  of  the  radome' 
ENDIF 


C      Prompt  for  Height  and  Radius  of  this  layer 
C 

WRITER,  *) 
WRITE(*,43) 
43         FORMAT(/8X,Please  enter  the  radome  layer  height,H,' 
1      /8X,'  and  the  outer  base  radius,R,  for  this  layer.') 
READ(*,*)  H(L),R(L) 

C      Prompt  for  shaping  function,  to  determine  equations  to  be  used 
C      to  generate  this  layer  of  the  radome 


C 


WRITE(V) 
WRITE(VO) 
10         FORMAT(/8X,'The  following  shaping  functions  are  available, ' 

1  /8X,'      1.  General  Ogive' 

2  /8X,'      2.  Tangent  Ogive' 

3  /8X,'      3.  von  Karman' 

4  /8X,'      4.  Power  Series' 

5  /8X,'      5.  Parabola,  curved' 

6  /8X,'      6.  Parabola,  capped  by  pointed  apex' 

7  //8X,'    Indicate  desired  shaping  function  by  entering  a  number' 

8  /10X,'  1,2,  ...or 6') 
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READ(*,*)  NSHAPE 

IF  (NSHAPE.EQ.1)  THEN 
C        Ogive 

SHAPE(L)  =  *Gen.  Ogive' 

CALL  OGIVE(L,X,Z,H,R,F,RESOLV,RESPTS  JMLMAX) 

ELSEIF  (NSHAPE.EQ.2)  THEN 
C        Tangent  Ogive 

SHAPE(L)  =  TAN  Ogive  ' 

CALL  TANOGIVE(L,X,Z,H,RJF,RESOLV  ,RESPTS,NLMAX) 

ELSEIF  (NSHAPE.EQ.3)  THEN 
C        von  Karman 

SHAPE(L)  =  'von  Karman' 

CALL  VONKRMN(L,XZ,H,R^ESOLVJlESPTSjvfLMAX) 

ELSEIF  (NSHAPE.EQ.4)  THEN 
C        Power  Series 

SHAPE(L)  =  Pwr  Series' 

CALL  POWER(L,X,Z,H,R  ,F,RESOLV,RESPTS  ,NLMAX) 

ELSEIF  (NSHAPE.EQ.5)  THEN 
C        Parabola  (curved) 

SHAPE(L)  =  Para  Curve' 

CALL  PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

ELSEIF  (NSHAPE.EQ.6)  THEN 
C        Parabola  with  apex 

SHAPE(L)  =  Para  APEX' 

CALL  APEX(L,XZ,HJ*F,RESOLVJ*ESPTS,NLMAX) 

ELSE 
C        Selection  out  of  limits 

WRITE  (*,*)  'Selection  is  not  available,  try  again.' 
WRITE  (*,*) 

114 


GOTO  39 
ENDIF 

C      Provide  listing  of  inner  layer  parameters 
C 

WRTTE(*,554) 

DO  11  I=1,L 
WRITE  (*,555)  I,SHAPE(I),H(I)^(I),(Faj)J=l,7) 

11  CONTINUE 

WRITE(*,*) 

12  CONTINUE 

C 

C  SHIFT  RADOME  SO  THAT  ORIGIN  IS  NEAR  THE  CENTER  OF  INNER  LAYER 

C    NOTE:  FACTOR  is  adjusted  for  optimal  placement  of  origin  for 

C        mesh  construction  by  EMRAD 

C 

FACTOR=0.4 

C      Max  inner  height  =  Z(NL,RESPTS) 
C 

SHFT=Z(  1  ,RESPTS)*(F  ACTOR) 

D0  88L=1,NL 
DO  77  1=1,  RESETS 
Z(L,I)=Z(LJ)-SHFT 
77  CONTINUE 

88         CONTINUE 


C 

Q  ****************************************************************** 

C     OUTPUT 

Q  ************************** ****************************** ********** 
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DO330L=l,NL 
C       OUTPUT  LAYER  X.Z  POINTS  FROM  TIP  TO  X-Z  PLANE 
DO220I=RESPTS,l,-l 
WRITE(2,*)  X(L,DZ(L,I) 
220         CONTINUE 

C      CLOSE  LAYER  BELOW  SHIFTED  X-Y  PLANE 

C       AT  Z=[-1*[SUM  OF  (LAYER  THICKNESS)] -SHFT] 

C 

C 

C        DETERMINE  LAYER  THICKNESS,  INNER  LAYER  HAS  NO  SHIFT  FOR  THICKNESS 

C         OTHER  LAYERS  HAVE  BOTTOM  THICKNESS  DETERMINED  BY  RADIAL  THICKNESS 

C         NOTE  -  THE  SHIFT  PERFORMED  FOR  THICKNESS  IS  ADDITIVE  FOR  LAYERS 

C 

IF(L.EQ.1)THEN 
THICK  =  0 

ELSE 
THICK  =  THICK  +  R(L)-R(L-1) 

ENDIF 
C 

BOTTM  =  (-1.0  *  THICK)  -  SHFT 

WRITE(2,*)  R(L),  BOTTM 

WRITE(2,*)  0.0,  BOTTM 

WRITE(2,*)  SEPRATN,SEPRATN 
330        CONTINUE 

C    SIGNAL  TERMINATION  OF  DATA  SET 
WRITE(2,100)  END,  END' 

C   Put  summary  of  data  to  screen 
C 

WRITE(*,554) 
DO440I=l,NL 
WRITE  (*,555)  I,SHAPE(I),H(I)^(I),(F(U)J=1,7) 
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440     CONTINUE 


C 

C      PROMPT  FOR  HEADER  INFORMATION 

C 

WRITE(V12) 

1 12  FORMAT(//7X,THE  HEADERS  ALLOW  THE  USER  TO  IDENTIFY  THIS7/7X, 
l'SET  OF  DATA  FROM  ALL  OTHER  SETS.', 

2//7X/PLEASE  ENTER  HEADER  #1  (64  CHARACTERS  MAX)7/7X) 

READ(VOO)  HDR1 

WRITE(M13) 

1 13  FORMAT(//7X,PLEASE  ENTER  HEADER  #2  (64  CHARACTERS  MAX)',//7X) 
READ(VOO)  HDR2 

C 

C      WRITE  HEADER  INFORMATION  INTO  MATERIAL  OUTPUT  DATA  FILE 

C 

WRITE(1,105)HDR1 

WRITE(1,105)HDR2 

C 

C      WRITE  LAYER  INFORMATION  INTO  MATERIAL  OUTPUT  DATA  FILE 

C 

WRITE(1,106)NL 

C     Prompt  for  material  Er  for  each  layer  of  construction 
C 

D0812L=1,NL 
WRITE(*,9)  L 

IF(L.EQ.1)THEN 

WRiTE(*,*) '       This  is  the  inner  core  of  the  radome' 
ELSEIF  (L.EQ.NL)  THEN 

WRITE(V) '       This  is  the  outer  layer  of  the  radome' 
ENDIF 
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C      Prompt  for  material  constants 
C 

WRITE(*,*)  "Enter  Er  for  this  Layer: ' 

READ(*,*)  ERR 

C     Ensure  that  Imag  Er  is  <  0 
C 

ERI  =  -0.000001 

ER(L)=CMPLX(ERR,ERI) 
UR(L)=CMPLX(1 .0,- 1  £-6) 
WRITE(1,102)  ER(L),UR(L) 
812  CONTINUE 


C       PLACE  COPY  OF  CONSTRUCTION  AT  BOTTOM  OF  MATERIAL  FILE 
C 

WRITE(  1,554) 
DO660I=l,NL 
WRITE  (1,555)  I,SHAPE(I),H(I)^(I),(F(U)J=1,7) 
660     CONTINUE 


:  —  FINAL  OUTPUT  TO  SCREEN 

WRITE(*,999)  MATFTL,DATFIL,DATFIL 
999  FORMAT(//8X, 
1'    Computations  and  output  are  completed.the  datafile  ,',A, 
2/.8X,'  holds  the  material  specifications  for  input  to  EMCAD.' 
3//,8X,'   The  datafile,  ',A,',  holds  the  radome  shape  data  for' 
4/,8X,'  viewing  in  the  CAD  package,  CURVE  DIGrTIZER. ' 
5//,8X,'   The  structure  file  required  for  input  to  EMCAD  can  be' 
6/.8X,'  generated  by  the  executing  the  program  EMCADIN,' 
7/,8X;  using,  ',A,' ,  as  input  to  EMCADIN, ' 
8//.8X,'   The  output  file  from  EMCADIN  will  be  in  proper  format' 
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9/.8X,'  for  input  to  EMCAD  for  the  structure  file.') 
WRITE(*,*) 
C 

lOOFORMAT(A) 

101  FORMAT(I5) 

102  FORMAT(4(E14.6)) 

103  FORMAT(I5,5(E14.6)) 

104  FORMAT(2(I5)) 

105  FORMATC  '.A) 

106  FORMAT(I5) 

554  FORMAT  (///////////////// 

1  //1X,LAYER\2X/  SHAPE    ,,2X/HEIGHT,2X,,RADIUS\ 

1  2X,'  Fl  *,2X,"  F2  '^X,'  F3  '^X,*  F4  ',2X,'  F5  ', 

2  2X,'  F6  ',2X,'  F7  71X,'- — '.2X,' ',2X,' ', 

3  2x,' — \2x; — ',2x,' — ',2x;-—'2x;-—\ 

4  2X,'-— ',2X,'— -\ix; — ') 

555  FORMAT  (1X,I5,2XA10,2X,F6.2,2X,F6.2J(F6.2)) 

9999  CONTINUE 

C 

C      END  PROGRAM  RADOME 

C 

CLOSE(l) 

CLOSE(2) 

STOP 

END 

*********************************** 

C  OGIVE  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  OGIVE(L,X,Z,H,R  ,F,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 

REAL  X(NLMAX,RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 

REALTHETA,C(5)^(5),H(5),DELZ(5),SHFT^ACTOR 

REAL  ALPHA1(5),  ALPHA(5),  B,  M,  RL.GAMMA,  F(NLMAX,7) 

COMPLEX  ER,UR 
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CHARACTER*  1  YN 

C    SET  CONSTANTS 
C 

PI=3.1415927 
DEGTORAD  =  PI/180. 
C 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 


C DEFINE  RADOME  LAYER  BY  OGIVE  EQUATIONS  .... 

DELZ(L)=H(L)/RESOLV 

C      Prompt  for  Shaping  factors  for  this  ogive  layer 
C 

WRITE(V) 
WRITE(*,43) 
43  FORMAT(8X,The  ogive  shaping  function  requires  three  shaping', 
l/8X,'factors  for  determining  the  radome  construction.' 
l//8X,'The  first  factor,  Fl,  determine  the  arc  of  curvature', 
l/8X,'of  the  radome  surface.  Please  enter  Fl:') 
READ(*,*)F(L,1) 

WRTTE(**) 
WRITE(*,45) 
45  FORMAT(8X,'  Two  other  shaping  factors  are  used  to  determine  the' 
l/8X,'relative  position  of  the  center  of  the  ogive  shaping', 
l/8X,'arc,  relative  to  the  radome  apex.', 
1//8X,'  The  first  of  these  factors,  F2,  determines  the  perp.  ', 
l/8X,'distance  from  the  z-axis  to  the  center  of  the  arc.', 
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1//8X,'  The  second  positioning  factor,  F3,  determines  the', 
l/8X,'the  distance  from  the  center  of  the  arc  to  the  apex', 
l/8X,'of  the  radome,  along  the  z-axis.', 
l//8X,TJlease  enter  the  position  factors  for  the  center  of  the', 
1/8X,'    curvature  arc,  F2  and  FS.*) 
READ(*,*)F(L,2),F(L,3) 

D0  33I=2,RESPTS-1 

Z(L,I)=Z(L,I- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=SQRT(  Fl**2-(F2-(H-z))**2 )  -  F3 

C 

X(L J)  =  SQRT(  F(L,1)**2  -( F(L,2)  -  (H(L)-Z(LJ))  )**2  )-F(L,3) 

33      CONTINUE 


C END  OF  RADOME  LAYER  .... 

RETURN 

END 

**********************£jvjj-)  Qp  ogive 

*********************************** 

C  TANGENT  OGIVE  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  TANOGrVE(L,XZ,H,R,F^ESOLV,RESPTS^LMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 

REAL  X(NLMAX,RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 

REAL  THETA,C(5),R(5),H(5),DELZ(5),SHFT,F  ACTOR 

REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL,GAMMA,  F(NLMAX,7) 

COMPLEX  ER,UR 

CHARACTER*  1YN 

C     SET  CONSTANTS 
C 

PI=3. 14 15927 

DEGTORAD  =  PI/180. 
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c 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C DEFINE  RADOME  LAYER  BY  TANGENT  OGIVE  EQUATIONS 

DELZ(L)=H(L)/RESOLV 

C      Determine  Shaping  factors  for  this  ogive  layer 
C 

F(L,4)  =  R(L)/2.0  +  2.0*(  H(L)**2.0  )/R(L) 

F(L,5)  =  F(L,4)  -  R(L)/4.0 

WRTTE(*,*) 

WRITE(*,43)F(L,4)JF(L,5) 
43  FORMAT(8X,'The  tangent  ogive  shaping  function  determines  two', 
l/8X,'factors  for  the  radome  construction.  These  are  calculated' 
l/SX.'based  on  the  height  and  radius  of  the  radome  layer.', 
l/8X,These  factors  are  F4  =  ',F6.2,'and  F5  =  \F6.2) 

D0  33I=2,RESPTS-1 

Z(L,I)=Z(L,I- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=SQRT(  F4**2  -  z**2  )  -  F5 

C 

X(L,I)  =  4.0*  (SQRT(  F(L,4)**2-  Z(L,I)**2  )  -  F(L,5)) 

33     CONTINUE 


C END  OF  RADOME  LAYER  .... 
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RETURN 
END 

****  ******************£j^jj-)  Qp  ffrj^  OGIVE 

C  von  Karman  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  VONKRMN(L,X,Z,H,R,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHA1  E,  NLMAX,  RESPTS,  RESOLV 

REAL  X(NLMAX,RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 

REALTHETA,C(5)J^(5)>H(5)X)ELZ(5),SHFT^ACTOR,SQPZF 

REAL  ALPHA1(5),  ALPHA(5),  B,  M,  RL.GAMMA 

COMPLEX  ER.UR 

CHARACTER*  1YN 

C    SET  CONSTANTS 
C 

PI=3. 14 15927 
SQP=SQRT(PI) 
DEGTORAD  =  PI/180. 
C 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 

C DEFINE  RADOME  LAYER  BY  VON  KARMAN  EQUATIONS  .... 

DELZ(L)=H(L)/RESOLV 

DO  33  I=2,RESPTS-1 

Z(L,I)=Z(L,M)+DELZ(L) 

C  EQUATION  OF  RADOME  IS 

C  x=(  R/(SQRT(PI)) )  *  SQRT(  FZ  -  SIN(2*FZ)/2  ) 
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:  for  FZ  =  ACOS(  1  -  2  *  (H-z)/H  ) 

FZ     =  ACOS(  1 .0  -  2.0  *  (  H(L)-Z(L,I) )  /  H(L) ) 
X(L  J)  *  R(L)/SQP  *  SQRT(  FZ  -  0.5*SIN(2*FZ)  ) 

33     CONTINUE 


C END  OF  RADOME  LAYER  .... 

RETURN 
END 

**********************END  op  VON  KARMAN 

*********************************** 

C  POWER  SERIES  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  POWER(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 

REAL  X(NLMAX,RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 

REALTHETA,C(5)^(5)^(5)X>ELZ(5),SHFTEACTOR 

REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL.GAMMA,  F(NLMAX,7) 

COMPLEX  ER,UR 

CHARACTER*  1YN 

C    SET  CONSTANTS 
C 

PI=3. 14 15927 
DEGTORAD  =  PI/180. 
C 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 


C DEFINE  RADOME  LAYER  BY  OGIVE  EQUATIONS 
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DELZ(L)=H(L)/RESOLV 

C      Prompt  for  Shaping  factor  for  this  power  series  layer 
C 

WRITE(*,*) 
WRITE(*,43) 
43  FORMAT(8X,'  The  power  series  shaping  function  requires  a', 
l/8X,'shaping  factor,  F6,  to  be  used  in  the  exponential  in', 
l/8X,'order  to  determine  the  radome  construction.', 
1//8X,'  When  F6  =  1.0  the  radome  shape  is  conical', 
1/8X,'  For  F6  =  0.75  the  radome  shape  approximates  the ', 
l/8X,T^ewtonian  contour  (minimizes  drag).  Please  enter  F6:') 
READ(*,*)  F(L,6) 

DO  33  I=2,RESPTS-1 

Z(L,I)=Z(L  J- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=R  *  (  (H-z)/H  )* *F6 

C 

X(LJ)  =  R(L)  *  ( (H(L)-Z(L,I))  /  H(L)  )**F(L,6) 

33     CONTINUE 


C END  OF  RADOME  LAYER  .... 

RETURN 
END 

**********************£j^|~)  Qp  POWER  SERIES 

*********************************** 

C  PARABOLIC  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS  JsfLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 

REAL  X(NLMAX,RESPTS),  Z(NLMAX RESPTS),  PI,  DELTHETA,  DEGTORAD 

REAL  THETA,C(5)^(5),H(5),DELZ(5),SHFTf  ACTOR 

REAL  ALPHA1(5),  ALPHA(5),  B,  M,  RL.GAMMA,  F(NLMAX,7) 
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COMPLEX  ER.UR 
CHARACTER*  1YN 

C    SET  CONSTANTS 
C 

PI=3.1415927 
DEGTORAD  =  PI/180. 
C 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 


C DEFINE  RADOME  LAYER  BY  PARABOLIC  EQUATION  .... 

C      Prompt  for  Shaping  factor  for  this  parabolic  layer 
C 

WPJTE(V) 
WRITE(*,43) 
43  FORMAT(8X,'  The  parabolic  shaping  function  requires  a', 
l/8X,'shaping  factor,  F7,  to  determine  the  curvature  of, 
l/8X,'the  parabolic  arc  for  the  radome  construction.' 
l//8X,'Please  enter  F7:') 
READ(*,*)F(L,7) 

C FOR  ROUNDED  PARABOLIC  RADOME  LAYER  .... 

C     ...  RADOME  IS  SCALED  TO  A  HEIGHT  =  SCALE  FUNCTION,C 

DELZ(L)=F(L,7)/RESOLV 

DO  33  I=2,RESPTS-1 
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Z(L,I)=Z(L,I- 1  )+DELZ(L) 

C  RECALL  EQUATION  OF  RADOME  IS  XA2=F7-Z,  X=SQRT(F7-Z) 

X(L,I)=SQRT(  F(L,7)-Z(LJ)  ) 


33     CONTINUE 

C 

C  SCALE  POINTS  TO  RADOME  DIMENSIONS 

C 

D0  44I=2,RESPTS-1 
X(L,I)=X(L,I)*R(L)/SQRT(  F(L,7) ) 
Z(L,I)=Z(L,I)*H(L)/F(L,7) 
44      CONTINUE 

C END  OF  ROUNDED  RADOME  LAYER  .... 


C END  OF  RADOME  LAYER  .... 

RETURN 
END 

**********************END  op  PARA  CURVE 

*********************************** 

C  POINTED  PARABOLIC  SHAPING  FUNCTION  SUBROUTINE 
SUBROUTINE  APEX(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 
INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 
REAL  X(NLMAX,RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA,C(5)3(5)JI(5),DELZ(5),SHFT,F  ACTOR 
REAL  ALPHA1(5),  ALPHA(5),  B,  M,  RL.GAMMA,  F(NLMAX,7) 
COMPLEX  ER,UR 
CHARACTER*  1YN 

C     SET  CONSTANTS 
C 
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PI=3. 1415927 
DEGTORAD  =  PI/180. 
C 

C    DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 
X(L,1)=R(L) 
Z(L,RESPTS)=H(L) 
X(L,RESPTS)=0.0 


C DEFINE  RADOME  LAYER  BY  POINTED  PARABOLIC  EQUATION  .... 


C      Prompt  for  Shaping  factor  for  this  parabolic  layer 
C 

WRITE(*,*) 
WRITE(*,43) 
43  FORMAT(8X,'  The  APEX  parabolic  shaping  function  requires  a' 
l/8X,'shaping  factor,  F7,  to  determine  the  curvature  of, 
l/8X,'the  parabolic  arc  for  the  radome  construction  and  the', 
1/8X, 'percent  of  the  radome  that  is  the  pointed  apex.', 
l//8X,'Please  enter  F7:') 
READ(*,*)F(L,7) 

C FOR  POINTED  RADOME  LAYER 


C     ...  RADOME  IS  SCALED  TO  A  HEIGHT  =  MODIFIED  SCALE  FUNCTION,F7+l 
DELZ(L)=(  F(L,7)+1.0  )/RESOLV 

DO  55  I=2,RESPTS-1 

Z(L,I)=Z(L,I- 1  )+DELZ(L) 
C 

IF(  Z(L,I).LE.(  F(L,7)  -1  ) )  THEN 
C  RECALL  EQUATION  OF  CONE  IS  XA2=F7-Z,  X=SQRT(F7-Z) 
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X(L,I)=SQRT(  F(L,7)-Z(L,I)  ) 
ELSE 
C  STRAIGHT  LINE  FORMULA  FOR  POINTED  END 

X(LJ)=(  (  F(L,7)  +1  )-Z(L,I)  )/2.0 
ENDIF 


55     CONTINUE 

C 

C  SCALE  POINTS  TO  RADOME  DIMENSIONS 

C 

D0  66I=2,RESPTS-1 
X(L,I)=X(L,I)*R(L)/SQRT(  F(L,7) ) 
Z(L,I)=Z(L,I)*H(L)/(  F(L,7)+L0  ) 
66    CONTINUE 

C END  OF  POINTED  RADOME  LAYER  .... 


C END  OF  RADOME  LAYER  .... 

RETURN 

END 

**********************£js^£)  Qp  para  APEX 
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Appendix  C 
CORTLD  PROGRAM  CODE 


PROGRAM  CORFLD 


C       Program  Created  by  LT  R  J.  Vince  (Aug-Dec  89) 
C       Modifications  by  M.A.  Morgan  (Nov-Dec  89) 

C       Computing  core  E-field  at  each  point  along  a  50x50  radome  grid 

C       in  the  radome's  field  coordinate  system. 

C 

INTEGER  BELL,FF,NPTS,L,J,K,K1  ,NA,NINC 

INTEGER  MSTART,MSTOP,N  1 ,NM 1  ,M,MF,MD,M  1  JVL2.MM 

INTEGER  I jrEJTM,KN,CA,CBJ4JNJSfLP,SELECT,SELCTl 

REAL  ALPHA(13)£>ALPHA(13),ANGINC(13) 

REAL  THETA,PHI,PI,DTR,RMIN,RCORE 

REALXP/iT,ZP,R,COST,SINT,COSP,SINP,S,C)CMP,SMP,SA,SAMX 

REALR1,R2^SQJEM,THJ1J«5TSLU1J>LOT(50,50),A(3,3) 

REAL  P(40),DP(40),TP(6)  J>P(6) 

REAL  EOTM,EOTE,EI(8),EDTERR(50,50,2),EDTER  1  .EDTER2 

COMPLEX  EXPltE(8)3DOT(50,50^)^DOTlJEDOT2 

COMPLEX  EP1  ,EP2,ET1  ,ET2,ER  1  £R2,JAY,RK,SRE,SRU,KR 

COMPLEX  ERAD(50,50,2),ETH(50,50,2),EPHI(50,50,2) 

COMPLEX  ER(  1 6),UR(  16),COEF(70,80),CJ(40),DCJ(40),CPLOT 

CHARACTER*  1  YN 

CHARACTER*8  PRTDAT 

CHARACTER*  12  PRTDT1,PRTDT2,PRTDT3 

CHARACTER*64  GRFLABJJDRl,HDR2TITLElTrTLE2,DUM 

$LARGE 

C 

C  Setting  constants 
BELL=CHAR(7) 
FF  =  CHAR(12) 
PI  =  3.14159265359 
DTR=  PI/180. 
JAY=  (0.0,1.0) 

C  Set  default  angle  of  rotation  of  planar  antenna  array  away  from  z  axis 

C  Set  default  aspect  ratio  for  plotting 
ASPECT  =1.35 

C  Set  number  of  points  of  resolution  across  array 
NPTS=50 
RPTS=(  FLOAT(NPTS)  -  1.  )/2. 

C**Initialize  arrays 

C  Initialize  the  E-field  matrices 
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c 

DO  1L=1,NPTS 
DO  1  J=1,NPTS 
DO  1  K=l,2 
ERAD(JJ-,K)=(0.,0.) 

ETH(J,L,K)=(0.,0.) 
EPHI(JJL,K)=(0.,0.) 
1    CONTINUE 

C  Prompt  for  input  data  filename  without  extension 
C 

WRITE(M4) 

READ(M01)  PRTDAT 

NCHAR=INDEX(PRTDAT,' ")  -1 

PRTDT1  =  PRTDAT(1:NCHAR)  //  '.INT 

OPEN(14JTLE=PRTDTl) 


C 


C 


TITLE  l=The  datafile  used  to  generate  fields  is ' 
TTTLE1=TTTLE1(1:40)//PRTDT1 

PRTDT2=  PRTDAT(1:NCHAR)  //  '.DAT' 
OPEN(  1 5  TTLE=PRTDT2) 


C  Write  header  from  structure  data  file 

C 

READ(14,101)HDR1 
READ{14,101)  HDR2 
READ(14,101)GRFLAB 
WRrrE(15,101)HDRl 
WRrrE(l  5,101  )HDR2 
WRITE(15,101)  GRFLAB 
WRITE(M01)HDR1 
WRITE(MOl)  HDR2 
WRITER,  10 1 )  GRFLAB 


C  Read  interior  core  Er,Ur 
C 

READ(14,507)  ER(1),DUM 

READH  4,507)  UR(1),DUM 

SRE=CSQRT(ER(1)) 

SRU=CSQRTfUR(l)) 

KR=SRE*SRU 

WRITEf*,*; 

wRrrEr*,*j'SRE=  ',sre 
wRrrEr*,*j'SRU=  ',sru 

WRITEfV/KR  =  ",KR 

WRITEriS ,*/  SRr^'.SRE,';  SRU=',SRU,';  KR  =  ',KR 

C  Read  minimum  diameter  of  missile  radome  core  expansion 
C      Assuming  in  Ko*r  Wavenumber  Units 

131 


c 

RE AD(  14,503)  RMTN.DUM 
RCORE=RMIN 

C  Read  number  of  incident  field  angles  (LE.  5) 
C 

READ(14,500)  NA.DUM 
NINC=2*NA 

C   Read  incident  angles 
C 

IPZ=0 

IMZ=0 

DOHI=l,NA 

READ(14,511)  DANGLE 

C       Checking  for  +/-  Z-Axis  Incident  Fields 
C 

IF(DANGLE.EQ.0.0)    IPZ=I 

IF(DANGLEEQ.  180.0)  IMZ=I 

DALPHA(I)=  DANGLE 

C      EMRAD  uses  Poynting  vector  angle  with  Z-axis  reference 
C       ANGle  of  incidence  =  180  -  DALPHA  from  EMRAD 
C 

ANGINC(I)=  180-  DANGLE 


WRTTE(*,*) '        Incoming  Inc.  Angle  #',I 

WRITE(*,*) '  EMRAD  Poynting  angle  of  incidence  =  ',DALPHA(I) 

WRITE(* ,*) '  Radome  incident  angle  =  ', ANGINC(I) 

WRITE(*,*) 

WRITE(15,*) '  Incoming  Inc.  Angle  #',I,'  =  \DALPHA(I) 

C  Set-up  array  of  angles  in  radian  measure 
C 

ALPHA(I)=DTR*DALPHA(I) 

1 1  CONTINUE 
NAI=1 
IF(NA.EQ.l)GOT013 

C 

C  Determine  angle  of  interest 

C 

12  CONTINUE 

WRTTE(*,*) '  Select  Angle  of  Interest  By  Number ' 

WRTTE(*,*) 

READ(*,*)  NAI 

IF  ((NAI.GE.l).OR.(NAI.LE.NA))  GO  TO  13 

WRITE(*,*) '  That  Number  is  Not  Available  ...  Try  Again' 

WRrTE(V) 
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GO  TO  12 
13  CONTINUE 

C 

C     Checking  if  Selected  Field  is  +/-  Z-Axis  Incident 
IM1=0 
IF((NAI.EQ.IPZ).OR.(NAI£Q.IMZ))  IM1=1 

WRITE(*,*) 

WRITE(*,*) '  The  angle  of  consideration  is  angle  #',NAI 

1  ,'  =  ',ANGINC(NAI),'  degrees* 

WRITER,*) 

WRrTE(15,*) '  The  angle  of  consideration  is  angle  #',NAI 

1  ,'  =  •,ANGINC(NAI),,  degrees' 

WRITE(15,*) 

C     Transformation  of  coordinates  to  gimble  mounted  planar  array 
C 

WRITE(*,*)'   The  planar  array  is  set  by  default  to  he  in  the' 
WRHfc(*,*)'    x-y  plane,  with  boresight  in  the  z-direction.' 
WRITE(*,*) 

WRITE(*,*)'    Do  you  want  to  change  the  array  boresight  ?' 
WRiTE(V)'  ( "Y"  or  "N" ) ' 

IFLG=0 
DTH=0.0 
READ(*,101)  YN 

rF(YN.EQ.T'.OR.YN.EQ.'y*)  THEN 
WRITE(*,*) 

WRTTE(*,*)  Enter  Boresight  Theta  Angle  (in  degrees)' 
READ(*,*)  DTH 
TH  =  DTH  *  DTR 
WRITE(**) 

WRITE(*,*)  Enter  Boresight  Phi  Angle  (in  degrees)' 
READ(V)  DFI 
H  =  DFI  *  DTR 

IFLG=1 

C    Loading  Transformation  Matrix  from  Array  to  Core  Coordinates 
C 

CT=COS(TH) 

ST=SIN(TH) 

CP=COS(FI) 

SP=siN(n) 

A(l,l)=CT*CP**2.O+SP**2.0 

A(1,2)=CP*SP*(CT-1.0) 

A(1,3)=ST*CP 

A(2,1)=CP*SP*(CT-1.0) 

A(2,2)=CT*SP**2.Oi-CP**2.0 

A(2,3)=ST*SP 

A(3,1)=-ST*CP 
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A(3,2)=-ST*SP 
A(3,3)=CT 

ENDIF 

C  Get  modal  information 

C 

READ(  14,505)  MSTART,MSTOP,DUM 

WRITEC*,*)'  MSTART=',MSTART,'  MSTOP=',MSTOP 

WRITE(15,*) "  MSTART=,,MSTART,'  MSTOP=',MSTOP 

C   Loop  modals  to  read  modal  coefficients 
C 

WRITEC*,*)  "Reading  Coefficients  and  Constructing' 

WRITEC*,*)  "Fields  for  Each  Azimuthal  ""m"  Mode' 

WRITEC*,*) 
C 

IF((MSTART.EQ.0).AND.(IM1.EQ.1))  THEN 

WRITEC*,*)  Terminating  M-Loop  at  m=l  for  +/-  Z-Axis  Incidence' 

WRITEC*,*) 

MSTOP=l 

ENDIF 

DO  777  MD=MSTART,MSTOP 

READ(14,500)  MJ)UM 

WPJTE(15,*)'M=',M 

FF(M.EQ.MD)  GO  TO  220 

WRITEC*,*)  "Error  in  Reading  "m"  -  Check  Data  File' 

STOP 

220  CONTINUE 
WRITEC*,*) 

WRITER,*)  'Azimuthal  Mode  "m"  =  ',M 

WRITEC*,*) 

IF(M.EQ.0)THEN 

EM  =1.0 

MF=1 
ELSE 

EM  =  2.0 

MF  =  0 
ENDIF 

WRITEC*,*)  "Reading  Coefficients' 

READ(14,500)  NM1.DUM 

D0  665ID=1,NA 

READ(14,500)  I.DUM 

IFCID.EQ.i)GOT0  221 

WRITEC*,*)  'Inc  Angle  Index  Read  Error  -  Check  Data  File' 

STOP 

221  CONTINUE 
WRrTE(15,*)'I=',I 
ITM=I 
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D0  222KN=lJSfMl 
CA=KN 
CB=CA+NM1 

READ{14,555)  COEF(CA,ITM),COEF(CB,ITM),DUM 
WRITE(15,*)  KN,COEF(CAJTM),COEF(CB  JTM),DUM 
222     CONTINUE 
ITE=I+NA 
D0  333KN=1,NM1 
CA=KN 
CB=CA+NM1 

READ(14,555)  COEF(CA,ITE),COEF(CB  JTE)X>UM 
WRITE(15,*)KN,COEF(CAjrE),COEF(CB,ITE),DUM 
333     CONTINUE 
665    CONTINUE 

C      Skipping  Field  Generation  for  m=0  if  +/-  Z-Axis  Incidence 
DF((M.EQ.0).AND.(TM1.EQ.1))  GO  TO  777 

WRITE(*,*)  'Generate  Fields  Across  Array' 

DX=RCORE/RPTS 
DY=DX 

D0  55J=1,NPTS 
D0  33L=lJSfPTS 

C        Determine  local  planar  array  coordinates 
C 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 

C       If  Array  Tilted  Then  Transform  Coordinates 

C       From  local  Planar  array  coordinates 

C  to  global  Core  coordinates.... 

C 

IF(IFLG£Q.1)THEN 
XC=A(1,1)*XP+A(1,2)*YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
ENDIF 

C        Determine  if  R  <=  RCORE,  Else  Skip  Field  Computation 
C 

RSQ=XC*XC+YC*YC+ZC*ZC 

R=SQRT(RSQ) 
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IF  (R.GT.RCORE)  GO  TO  32 

C         Determine  spherical  Bessel  function  scaled  radius  argument 
C 

RK=KR*R 

C  Determine  Complex  Spherical  Bessel  function  and  Derivative 

C 

NJN=M+NM1+MF 

CALL  CSBSL(RK,NJN,CJ  ,DCJ) 

C      Compute  Core  Coordinates  and  Functions 
C 

THETA=ACOS(ZC/R) 

PHI=ARCTAN2(YC,XQ 

C=COS(THETA) 

S=SIN(THETA) 

CMP=EM*COS  (M*PHI) 

SMP=EM*SIN(M*PHI) 

C       Determine  Legendre  Polynomial  and  Derivative 
C 

NLP=NM1+MF 

CALL  LPAD(MJSfLP,C,PJ)P) 

C      Consider  only  the  incident  angle  of  interest 
C 

I=NAI 

ITM=I 

ITE=ITM+NA 

C      Initialize  additive  KN-sum  fields  to  zero 
C 

EP1=(0.,0.) 

EP2=(0.,0.) 

ET1=(0.,0.) 

ET2=(0.,0.) 

ER1=(0.,0.) 

ER2=(0.,0.) 


D0  22KN=1,NM1 
N=KN+M+MF-1 
NJ=N+1 
NL=KN+MF 
CA=KN 
CB=CA+NM1 

ET1=ET1+C0EF(CA,ITM)*M*CJ(NJ)*P(NL) 
1  +COEF(CB,ITM)*DCJ(NJ)*DP(NL) 

EP1=EP1+C0EF(CA,ITM)*CJ(NJ)*DP(NL) 
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1  +COEF(CBXTM)*M*DCJ(NJ)*P(NL) 

ERl=ERl+COEF(CBXTM)*N*(N+l)*CJ(NJ)*P(NL) 

ET2=ET2+COEF(CA,ITE)*M*a(NJ)*P(NL) 
1  +COEF(CBXTE)*DCJ(NJ)*DP(NL) 

EP2=EP2+COEF(CAJTE)*CJ(NJ)*DP(NL) 
1  +COEF(CBJTE)*M*DCJ(NJ)*P(NL) 

ER2=ER2+COEF(CBJTE)*N*(N+l)*a(NJ)*P(NL) 

22         CONTINUE 


ETH(J  ,L,  1  )=ETH(J,L,  1  >J  AY*SRU*CMP*ET1/R 

EPHI(J,L,1)=EPHI(J,L,1)+JAY*SRU*SMP*EP1/R 

ERAD(J^,1)=ERAD(J,L,1)-JAY*S*CMP*ER1/(SRE*RSQ) 

ETH(J^^)=ETH(J,L,2>+SRU*SMP*ET2/R 

EPHI(J,L,2)=EPHI(JtL,2)+SRU*CMP*EP2/R 

ERAD(J,L,2)=ERAD(J,L,2)+S*SMP*ER2/(SRE*RSQ) 


32  CONTINUE 

33  CONTINUE 

* 

55   CONTINUE 
C 
C=  End  of  loops  across  antenna  plane : 


777  CONTINUE 
C 
C**  End  of  modal  loop  for  field  analysis 


>*************** 


C    Determine  magnitude  and  phase  of  Computed  field  in  the  direction 
C    of  the  theoretical  field  vector,  and  the  magnitude  of  the  error. 
C 

SINA=SIN(ALPHA(NAI)) 
COSA=COS(ALPHA(NAI)) 


C   Loop  through  array  to  determine  field  components 
C 

D0  634J=1,NPTS 

D0  633L=1,NPTS 

C        Determine  local  planar  array  coordinates 
C 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 
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C       If  Array  Tilted  Then  Transform  Coordinates 

C       From  Array  to  Core .... 

C 

IF(IFLG£Q.1)THEN 
XC=A(1,1)*XP+A(1,2)*YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 

ELSE 

XC=XP 

YOYP 

ZC=ZP 
ENDIF 

RSQ=XC*XC+YC*YC+ZC*ZC 

R=SQRT(RSQ) 

C        Determine  if  R  <=  RCORE,  Else  Zero  Fields 
C 

IF  (R.LE.RCORE)  THEN 

C      Compute  Core  Coordinates  and  Functions 
C 

THETA=ACOS(ZC/R) 

PHI=ARCTAN2(YC,XC) 

CT=COS(THETA) 

ST=SIN(THETA) 

CP=COS(PHI) 

SP=SIN(PHI) 

C      NOTE:  R  IS  GIVEN  IN  TERMS  OF  KO*R 

C  ==  KO  IS  INCLUDED  IN  R,  DEFINE  KO=l .0 

C 


C 


KO=1.0 


C      TM  and  TE  Plane  waves  are  of  unit  magnitude 
C 

EOTE=1.0 

EOTM=1.0 


C    Complex  incident  field  is  exponential(-j*ko*p.r) 

C     use  field  direction  for  incident  angle  180  to  90  degrees 


EXP1=  CEXP(  -JAY*KO*(  XC*SINA+ZC*COSA  ) ) 


C      Compute  plane  wave  field  unit  vectors  and  components 
C 

C       E  TM  =  cos(alpha)Ux  -  sin(alpha)Uz 
C 

C        ETMT 
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EI(1)  =  COSA*CT*CP  -  SINA*(-ST) 

E(1)  =  E0TM*EXP1*EI(1) 
C        ETMP 

EI(2)  =  COSA*(-SP) 

E(2)  =  E0TM*EXP1*EI(2) 
C        ETMR 

EI(3)  =  COSA*ST*CP  -  SINA*(CT) 

E(3)  =  E0TM*EXP1*EI(3) 
C      E  TM  TOT 

E(4)=SQRT(CABS(E(1))**2.0+CABS(E(2))**2.0+CABS(E(3))**2.0) 

C     E  TE  =  Uy 
C 

C        ETET 

EI(5)  =  CT*SP 

E(5)  =  EOTE*EXPl*EI(5) 
C        ETEP 

EI(6)  =  CP 

E(6)  =  EOTE*EXPl*EI(6) 
C        ETER 

EI(7)  =  ST*SP 

E(7)  =  EOTE*EXPl*EI(7) 
C        E  TE  TOT 

E(8)=SQRT(CABS(E(5))**2.0+CABS(E(6))**2.0+CABS(E(7))**2.0) 

********** 

C   Find  magnitude  and  phase  of  Computed  field  in  the  direction 

C    of  the  theoretical  field  vector 

C      E.ei/Ei=  ( eth*ETH+ephi*EPHI+erad*ERAD  )/EO*EXPl 

C 

ED0T1=EI(1)*ETH(J^,1)+EI(2)*EPHI(J^,1)+EI(3)*ERAD(J^,1) 

EDOT(J^,l)=EDOTl/(EOTM*EXPl) 

EIX)T2=EI(5)*ETH(JJL,2)+EI(6)*EPffl(JJ.,2)+EI(7)*ERAD(J^,2) 
EDOT(JJL,2)=EDOT2/(EOTE*EXPl) 

C   Find  rel.  magnitude  of  Computed  field  perpendicular  to  the  direction 

C    of  the  theoretical  field  vector,  the  error  component  in  the 

C    Computed  field. 

C 

C    Rel  error  =  I E  -  (E.ei)  ei  I  /  EO 

C 

C    =  sqrt{  Eth  -  (E.ei)*ethl**2  +  Erad  -  (E.ei)*eradl**2 

C  +  Ephi  -  (E.ei)*ephil**2  }/  EO 

C 

C   Field  error  magnitude  is  computed 

C 

EDTER1=  CABS(ETH(J,L,1)  -EDOT1*EI(1))**2.0 

1       +CABS(EPffl(J,L,l)-EDOTl*EI(2))**2.0 

1       +  CABS(ERAD(J,L,1)-EDOT1*EI(3))**2.0 
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EDTERR(J,L,1)=SQRT(EDTER1)/(E0TM) 

EDTER2=  CABS(ETH(J,L,2)  -EDOT2*EI(5))**2.0 
1       +  CABS(EPHI(J,L,2)-EDOT2*EI(6))**2.0 
1       +  CABS(ERAD(J,L,2)-EDOT2*EI(7))**2.0 
EDTERR(JJL,2)=SQRT(EDTER2)/(EOTE) 


ELSE 

C         R  >  RCORE  ,  Set  field  components  to  zero 
C 

EDOT(J,L,1)=0.0 

EDOT(J,L,2)=0.0 

EDTERR(J,L,1)=0.0 

EDTERR(J,L,2)=0.0 
ENDIF 

633  CONTINUE 

634  CONTINUE 

WRTTE(*,*)  'SPHERICAL  FIELDS  COMPUTED' 
WRTTE(*,*) 
C   Optional  Listing  of  Coordinates  and  E-Fields 
C 

WRTTE(V)  Check  Values  of  Coordinates  and  E-Field  ?  (Y/N): 
READ(MOl)  YN 

IF((YN.EQ.T').OR.(YN.EQ.y))  THEN 
757     CONTINUE 

WRITE(*,*)  "Enter  Array  Triple  Index  to  Check  (J,L,I)' 

WRITE(*,*) '        J=  x-point  1,2,..,\NPTS 

WRITE(*,*) '        L=  y-point  1,2,..,',NPTS 

WRITE(*,*) '        1=  1  for  TM,  2  for  TE' 

WRITE(*,*) 

READ(*,*)J,L,I 

EPM=C  ABS(EPHI(J  L,I)) 

ETM=CABS(ETH(J,L,I)) 

ERM=CABS(ERAD(J,L,I)) 

ETOT=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

C        Determine  local  planar  array  coordinates 
C 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 

C       If  Array  Tilted  Then  Transform  Coordinates 
C       From  Array  to  Core  .... 
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c 

IF(IFLG£Q.1)THEN 

XC= A(  1 , 1  )*XP+A(  1 ,2)*  YP+ A(  1 ,3)*ZP 

YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 

ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 
ELSE 

XC=XP 

YC=YP 

ZC=ZP 
ENDIF 

WRITE(*,*)  *XP,YPZP:',XP,YPZP 
WRITE(*,*)  ,XC,YC,ZC:,,XC,YC,ZC 
WRITE(V) 

WRITE(*,*)  'IE-RADI  =',ERM 
WRITE(*,*)  'E-THI  =',ETM 
WRITE(*,*)  "E-PHII  ='EPM 
WRITEC*,*)  1E-T0TI  ="£TOT 
WRITE(*,*) 

WRITE(*,*)  "Look  at  Another  Point  ?  (Y/N):' 
READ0UO1)  YN 

IF((YN.EQ.  Y*).OR.(YN.EQ.y))  go  to  757 
ENDIF 

Q ********************************************************************** 

C     Field  representation  selection  menu 

Q  *********************************** 

C 

700  WRITER  ,701) 

Q  ****************************************************************** 

701  format(////////////5X,'CORFLD  OUTPUT  GRAPHICS  PRESENTATION^X, 
l'CORE  FIELD  REPRESENTATION  SELECTION  MENU7/5X, 

2Please  select  the  type  of  field  representation  of  interest 

3V5X, 

4'CORE  will  display  the  following  field  representauons:7/8X, 

5'1.    TM  INCIDENT  FIELD,    E-THETA78X, 

6'2.    TM  INCIDENT  FIELD,    E-PHI  78X, 

7'3.    TM  INCIDENT  FIELD,    E-RADIAL  'JSX, 

8'4.    TM  INCIDENT  FIELD,    E-TOTAL  FIELD  MAGNITUDE  7/8X, 

9'5.    TE  INCIDENT  FIELD,    E-THETA78X, 

1'6.    TE  INCIDENT  FIELD,    E-PHI  'JSX, 

27.    TE  INCIDENT  FIELD,     E-RADIAL  78X, 

3'8.    TE  INCIDENT  FIELD,    E-TOTAL  FIELD  MAGNITUDE  7/8X, 

4'9.    COMPARE  COMPUTED  FIELD  TO  INCIDENT  FTELD7/8X, 

5' 10.    CHANGE  ASPECT  RATI07/8X, 

6'0.    FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE7/5X, 

6'Indicate  your  selection  entering  "l","2",...,"10",or  "0" 

7',5X) 

Q  ********************************** 

READ(*,*)  SELECT 
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C     Check  response  to  menu 
IF(SELECT£Q.1)THEN 

Q  ****************************************************************** 

C        DISPLAY  TM  INCIDENT  FIELD,  E-THETA 

C       CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,715) 

Q  ****************************************************************** 

715  format(////////////5X,*CORE  FIELD  PRESENTATIONS75X, 
ITou  have  selected  the  following  field  element' ,//5X, 
21.    TM  INCIDENT  FIELD,    E-THETA7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNnTJDE,7/5X, 
4Please  select  the  type  of  field  representation  of  interest 
5'J8X, 

61.    PHASE'^SX, 
7'2.    MAGNTrUDE78X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
97/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C     ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  712J=1,NPTS 
D0  711L=1,NPTS 
CPLOT=ETH(J,L,l) 
IF  (SELCT1.EQ.1)  THEN 

PLOT(JX)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 

PLOT(J,L)=CABS(CPLOT) 
ENDIF 

711  CONTINUE 

712  CONTINUE 
C 

IF(SELCT1.EQ.1)THEN 

TITLE2  =  ' TM  INCIDENT  FIELD,    E-THETA   PHASE' 
ELSE 

TTTLE2  =  '  TM  INCIDENT  FIELD,     E-THETA    MAGNITUDE 
ENDIF 
C 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.2)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TM  INCIDENT  FIELD,  E-PHI 

C       CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,725) 

Q  ****************************************************************** 
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725  format(////////////5X,'CORE  FIELD  PRESENTATIONS',^, 
l'You  have  selected  the  following  field  element7/5X, 
22.    TM  INCIDENT  FIELD,    E-PHI7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNrrUDE,7/5X, 
4  Please  select  the  type  of  field  representation  of  interest 
5-./8X, 

6'1.   PHASE'^SX, 
7'2.    MAGNrrUDE'y8X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
97/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C     ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  722J=1,NPTS 
D0  721L=1,NPTS 
CPLOT=EPHI(J,L,l) 
IF(SELCT1£Q.1)THEN 

PLOT(J0-)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 

721  CONTINUE 

722  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TTTLE2  =  '  TM  INCIDENT  FIELD,    E-PHI     PHASE    ' 
ELSE 

TTTLE2  =  ' TM  INCIDENT  FIELD,    E-PHI     MAGNITUDE" 
ENDIF 


Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.3)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TM  INCIDENT  FIELD,  E-RADIAL 

C        CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,735) 

Q  ****************************************************************** 

735  format(////////////5X,'CORE  FIELD  PRESENTATIONS75X, 
l'You  have  selected  the  following  field  element' y/5X, 
2'3.    TM  INCIDENT  FIELD,     E-RADIAL7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,7/5X, 
4Please  select  the  type  of  field  representation  of  interest 
5',/8X, 

6'1.   PHASE'^SX, 
7'2.    MAGNITUDE78X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
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97/5X) 
p     ********************************** 

REAEK*,*)  SELCT1 

C     ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  732J=1,NPTS 
D0  731L=1,NPTS 
CPLOT=ERAD(J,L,l) 
IF  (SELCT1.EQ.1)  THEN 

PLOT(JtL)=(ARCTAN2(IMAG(CPLOT)^EAL(CPLOT)))/DTR 
ELSE 

PLOT(J,L)=CABS(CPLOT) 
ENDIF 

731  CONTINUE 

732  CONTINUE 
C 

IF(SELCT1.EQ.1)THEN 

TITLE2  =  '  TM  INCIDENT  FIELD,     E-RADIAL  PHASE  ' 
ELSE 

TITLE2  =  '  TM  INCIDENT  FIELD,     E-RADIAL  MAGNITUDE' 
ENDIF 
C 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.4)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 
C      ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  742J=1,NPTS 

D0  741L=1,NPTS 

EPM=CABS(EPHI(J  JL,  1 )) 
ETM=CABS(ETH(J  i,  1 )) 
ERM=CABS(ERAD(J,L,1)) 
PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

741  CONTINUE 

742  CONTINUE 

TITLE2  =  ' TM  INCIDENT  FIELD,     E-TOTAL  FIELD  MAGNITUDE  ' 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.5)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TE  INCIDENT  FIELD,  E-THETA 

C        CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,755) 

r  ****************************************************************** 
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755  format(////////////5X,'CORE  FIELD  PRESENTATI0NSV5X, 
lTou  have  selected  the  following  field  element' ,//5X, 
2'5.    TE  INCIDENT  FIELD,    E-THETA7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNTTUDE,7/5X, 
4  Please  select  the  type  of  field  representation  of  interest 
578X, 

6*1.    PHASEV8X, 
77.    MAGNITUDE'ySX, 
81ndicate  your  selection  entering  "  1 "  or  "2" 
97/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C     ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  752J=1,NPTS 
D0  751L=1,NPTS 
CPL0T=ETH(JJL2) 
IF(SELCT1.EQ.1)THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 

PLOT(J,L)=CABS(CPLOT) 
ENDIF 

751  CONTINUE 

752  CONTINUE 
C 

IF(SELCT1£Q.1)THEN 

TTTLE2  =  *  TE  INCIDENT  FIELD,    E-THETA   PHASE' 
ELSE 

TTTLE2  =  ' TE  INCIDENT  FTELD,     E-THETA    MAGNITUDE' 
ENDIF 


p        ****************************************************************** 
ELSEIF  (SELECT.EQ.6)  THEN 

£  ****************************************************************** 

C       DISPLAY  TE  INCIDENT  FTELD,  E-PHI 

C       CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,765) 

Q  ****************************************************************** 

765  format(////////////5X,'CORE  FIELD  PRESENTATIONS75X, 
lTou  have  selected  the  following  field  element7/5X, 
2'6.    TE  INCIDENT  FTELD,    E-PHI7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,7/5X, 
4Please  select  the  type  of  field  representation  of  interest 
578X, 

6'1.    PHASE'^SX, 
7'2.    MAGNITUDE78X, 
8'Indicate  your  selection  entering  "  1 "  or  "2" 
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97/5X) 

Q  ********************************** 

READ(V)  SELCT1 

C      ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  762J=1,NPTS 
D0  761L=1,NPTS 
CPLOT=EPHI(JtL,2) 
IF  (SELCT1.EQ.1)  THEN 

PLOT(JLMARCTAN2(TMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 

PLOT(JJL)=CABS(CPLOT) 
ENDIF 

761  CONTINUE 

762  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TrTLE2  =  '  TE  INCIDENT  FIELD,     E-PHI      PHASE  ' 
ELSE 

TITLE2  =  ' TE  INCIDENT  FIELD,     E-PHI      MAGNITUDE  ' 
ENDIF 
C 


C 


****************************************************************** 


ELSEIF  (SELECT.EQ.7)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TE  INCIDENT  FIELD,  E-RADIAL 

C        CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,775) 

Q  ****************************************************************** 

775  format(////////////5X,'CORE  FIELD  PRESENTATIONS75X, 
l'You  have  selected  the  following  field  element' ,//5X, 
27.    TE  INCIDENT  FIELD,     E-RADIAL7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,7/5X, 
4Please  select  the  type  of  field  representation  of  interest 
578X, 

6'1.    PHASE'^X, 
7'2.    MAGNITUDE78X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
97/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C      ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  772J=1,NPTS 
D0  771L=1,NPTS 
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CPL0T=ERAD(JJL,2) 
IF  (SELCT1.EQ.1)  THEN 

PL0T(J,L)=(ARCTAN2(IMAG(CPL0T),REAL(CPL0T)))/DTR 
ELSE 

PLOT(JX)=CABS(CPLOT) 
ENDIF 

771  CONTINUE 

772  CONTINUE 
C 

IF(SELCT1£Q.1)THEN 

TITLE2  =  '  TE  INCIDENT  FIELD,     E-RADIAL  PHASE  ' 
ELSE 

TTTLE2  =  '  TE  INCIDENT  FIELD,    E-RADIAL  MAGNTrUDE' 
ENDIF 
C 

Q  ***********************•***=(<********************  +  ************  +  **»  +  * 

ELSEIF  (SELECT.EQ.8)  THEN 

Q  ****************************************************************** 

C        DISPLAY  TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

D0  782J=1,NPTS 

D0  781L=1JSTPTS 

EPM=CABS(EPHI(Ji,2)) 

ETM=CABS(ETH(J,L,2)) 

ERM=CABS(ERAD(J,L,2)) 

PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

781  CONTINUE 

782  CONTINUE 

TITLE2  =  ' TE  INCIDENT  FIELD,    E-TOTAL  FIELD  MAGNITUDE  ' 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.9)  THEN 

Q  ********************************************************************** 

c 

C     Incident  Field  comparison  representation  selection  sub-menu 

Q      ************************************************************ 

C 

800  WRITE(*,801) 

Q  ****************************************************************** 

801  format(////////////5X,'CORE  COMPARISON  TO  IDEAL  FTELD7/5X, 
l'CORE  FIELD  REPRESENTATION  SELECTION  MENU  27/5X, 
2Please  select  the  type  of  field  representation  of  interest 

375X, 

4'CORE  will  display  the  following  field  representations:7/8X, 

5*1 .    TM  INCIDENT  FIELD,    Computed  Field  dot  Ideal  Field78X, 

6'2.    TM  INCIDENT  FIELD,    Scaled  Error  component  7//8X, 

73.    TE  INCIDENT  FIELD,    Computed  Field  dot  Ideal  Field',/8X, 

8'4.    TE  INCIDENT  FIELD,    Scaled  Error  component  7//8X, 

9'Indicate  your  selection  entering  "1","2","3",  or  "4" 
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17/5X) 

Q  ********************************** 

READ(*,*)SELCT2 

C     Check  response  to  menu 

****** 

IF  (SELCT2.EQ.1)  THEN 

Q  ****************************************************************** 

C       DISPLAY  TM  INCIDENT  FIELD,    Computed  Field  dot  Ideal  Field 

C       CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 

WRITE(*,815) 

p      ****************************************************************** 

815  format(////////////5X,'CORE  FIELD  PRESENTATIONS',^, 
lTou  have  selected  the  following  field  element'  J/5X, 
2'1 .   TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE, V/5X, 
4Please  select  the  type  of  field  representation  of  interest 
578X, 

6*1.   PHASE',/8X, 
7'2.    MAGNTTUDEV8X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
9V/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C      ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  812J=1,NPTS 
D0  811L=1,NPTS 
CPLOT=EDOT(J,L,l) 
IF(SELCT1£Q.1)THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 
PLOT(J,L)=CABS(CPLOT) 
ENDIF 

811  CONTINUE 

812  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TTTLE2=TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE' 
ELSE 

TITLE2=TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN.' 
ENDIF 
C 


C 


****************************************************************** 

ELSEIF  (SELCT2.EQ.2)  THEN 

****************************************************************** 
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C       DISPLAY  TM  INCIDENT  FIELD,    Scaled  Error  component 

D0  822J=1,NPTS 
D0  821L=lJSfPTS 
PLOT(J,L)=EDTERR(J,L,l) 

821  CONTINUE 

822  CONTINUE 

TTTLE2  =  TM  INCIDENT  FIELD,    Scaled  Error  component  ' 

Q  ****************************************************************** 

ELSEIF  (SELCT2.EQ.3)  THEN 

Q  **************************************************************** 

C       DISPLAY  TE  INCIDENT  FIELD,    Computed  Field  dot  Ideal  Field 
C       CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WPJTE(*,815) 

Q  ****************************************************************** 

835  format(////////////5X,'CORE  FIELD  PRESENTATIONS',^, 
lYou  have  selected  the  following  field  element7/5X, 
2'3.   TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field7/8X, 
3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE.7/5X, 
4  "Please  select  the  type  of  field  representation  of  interest 
578X, 

6'1.    PHASE' JSX, 
7*2.    MAGNITUDE78X, 
8'Indicate  your  selection  entering  "1"  or  "2" 
97/5X) 

Q  ********************************** 

READ(*,*)  SELCT1 

C      ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

D0  832J=1,NPTS 
D0  831L=1,NPTS 
CPLOT=EDOT(J,L,2) 
IF(SELCT1.EQ.1)THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 
ELSE 

PLOT(J,L)=CABS(CPLOT) 
ENDIF 

831  CONTINUE 

832  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TITLE2=TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE' 
ELSE 

TITLE2=TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN.' 
ENDIF 
C 

p        ************************* ***************************************** 
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ELSEIF  (SELCT2.EQ.4)  THEN 

Q  ****************************************************************** 

C       DISPLAY  TE  INCIDENT  FIELD,    Scaled  Error  component 

D0  842J=1,NPTS 
D0  841L=lJsrPTS 
PLOT(J,L)=EDTERR(J,L,2) 

841  CONTINUE 

842  CONTINUE 

TTTLE2  =  TE  INCIDENT  FIELD,    Scaled  Error  component  ' 

Q  ****************************************************************** 

ELSE 

Q  ****************************************************************** 

C       ( (SELCT2JLT.0).OR.(SELCT2.GT.4) ) 
GO  TO  800 
ENDIF 
C  — 

Q ********************************************************************** 


ELSEIF  (SELECT.EQ.10)  THEN 

p        ****************************************************************** 

C       ASPECT  RATIO 

791  write(*,792) 

792  FORMAT(/5X;**********  ASPECT  RATIO  **********7/5X, 

1'IN  ORDER  TO  AVOID  DISTORTION  ON  THE  SCREEN,  IT  IS  NECESSARY75X, 
2TO  DEFINE  THE  NUMBER  OF  ROWS  PER  INCH  AND  THE  NUMBER  OF75X, 
3'COLUMNS  PER  INCH.  THESE  DEFINITIONS  CHANGE  AS  THE75X, 
4"RESOLUTION  OF  THE  SCREEN,  AS  WELL  AS  THE  TYPE  OF  SCREEN75X, 
5'CHANGE.7/5X, 
C      'THIS  PROGRAM  WAS  WRITTEN  IN  MODE  16  (EGA)  USING75X, 
C     6'A  GB-1  VIDEO  BOARD  AND  A  NEC  MULTISYNC  MONITOR.  THE75X, 
C     7'OPTIMUM  ASPECT  RATIO  FOR  THIS  CONFIGURATION  WAS  .657//5X, 
8PLEASE  INPUT  ASPECT  RATI0.7/5X) 

read(*,*)  aspect 
goto  700 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.O)  THEN 

Q  ****************************************************************** 

C       FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 
goto  999 

ELSE 

Q  ****************************************************************** 

C       ( (SELECT.LT.0).OR.(SELECT.GT.9) ) 
GO  TO  700 
ENDIF 
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Q     ******+****************++**+++*++*+**++++*++++** +*++******************* 


WRITE  (*,*) '  CALL  PLOTTING  ROUTINE" 
CALL  PLOT3D(PLOT,50t50ASPECT,GRFLAB, 
+  TITLE  l,TTTLE2,HDRlJIDR2tANGINC(NAI),DTH) 

GO  TO  700 

C  End  of  graphic  output  routine 

999  CONTINUE 

C 

C      Program  completed 
C 

C    I/O  FORMAT  STATEMENTS 
C 
10  FORMAT(/////////////////7X/*****  WELCOME  TO  CORFLD  ******  J/7X, 
1THIS  PROGRAM  COMPUTES  INTERIOR  EM  FIELD  WITHIN  PENETRABLE'/7X, 
2'BODIES  OF  REVOLUTION  USING  THE  FIELD  COMPONENT  rNFORMATION'/7X, 
3'OUTPUT  BY  THE  PROGRAM  EMCAD//7X, 
4PLEASE  PRESS  ANY  KEY  TO  CONTINUE.',/) 
14  FORMAT(//////////////////////7X, 
1THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT/7X, 
2PESULTS  OF  ALL  NUMERICAL  CALCULATIONS  CONDUCTED  BY  EMCAD.',/7X, 
3THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO  EXTENSION  IS*/7X, 
4'REQUIRED  OR  DESIRED.  EMCAD  AUTOMATICALLY  APPENDS  AN  'J7X, 
SPXTENSION  OF  .INT  TO  YOUR  FILENAME.7/7X, 
6PLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILE.7/7X) 
23  FORMAT(/7X,ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  \UJPX) 
25  FORMAT(//////7X,PLEASE  PRESS  ANY  KEY  TO  CONTINUE. 7///) 

100  FORMAT(/  *********************  CORE  OUTPUT  DATA  **************** 

101  FORMAT(A) 

102  FORMAT(I5) 

103  FORMAT(4(E14.6)) 

104  FORMATC  ***  PROGRAM  ABORTED  BECAUSE ***') 

107  FORMAT(//7X,I5,2(2X,'(\lPEl  1.3^X,1PE1 1.3,')')) 

109  FORMAT(//7X,TNCIDENT  FIELD  ANGLES7/9X,,N,,10X,■ALPHA(N)•) 
110FORMAT(//7X,I3,7XJ;5.0,,  DEC) 

111  FORMAT(/7X,'COMPLEX  Er(n)  AND  Ur(n)") 

1 12  FORMAT(//7X,PNTER  INTERIOR  DATA  FILE  (DrFTLENAME.EXTENSION):  ') 
1 15  FORMAT(//7X,TiNTER  CAPTION  OR  LABEL     ') 

124  FORMAT(//7X,TNDEXING  PROGRAM  THROUGH  VALUES  OF  M') 

125  FORMAT(7X,'M-LOOP ....  M  =  \I5) 

135  FORMAT(7X,PX  M-LOOP) 

136  FORMAT(//7X,'*********  CORE  ANALYSIS  COMPLETED  **************■) 
204  FORMAT(//7X  '  **************  M  =  '  15  '  **************va 

300  FORMAT(7X,I5,5E12.3) 
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500  FORMAT(3X,I5,A64) 

503  FORMAT(3X,E14.6,A64) 

505  FORMAT(3X,I5,2XJ5A64) 

507  FORMAT(3X,lPEl  1.3.2X.1PE1 1.3.2X.A64) 

510  FORMAT(3X,F5.0,A64) 

511FORMAT(3XJF5.0) 

555FORMAT(3X,E12.63X,E12.6,4X,E12.6,3X,E12.6,3X,A64) 

566  FORMAT(3X,A10,E12.6,3X,E12.6) 

567  FORMAT(3X,A  10.E  12.6) 

576  FORMAT(3X,A5,I5,A2,E12.6,3X,E12.6) 

577  FORMAT(3X,A5312.6) 

578  FORMAT(3X,E12.6) 


99  CONTINUE 
C 

C    Program  termination 
C 

WRITE(*,136) 

STOP 

END 

**^*^*******^r******** ******************************  *********** 

SUBROUTINE  CSBSL  (Z,N1,CJ,DCJ) 

COMPLEX  Z,Z2,F1,F2,FST0,ALPHA)CJ(*),DCJ(*) 
C 

C     This  routine  computes  Riccati  spherical  bessel  functions  and 
C     derivatives  for  complex  Z  using  downward  recurrence  as  in 
C     Miller's  algorithm,  but  with  highly  accurate  starting  values 
C     using  truncated  Tatlor's  series. 
C 

C     Returns  CJ(n)=J[n-l]  (kr)  for  n=l..JSfl 
C  DCJ(n)=J'[n-l]  (kr) 

C 

C     DEFINING  CONSTANTS 
C 

ZM2=(CABS(Z))**2 

NMX=N1+ZM2 

Z2=(Z*Z)/2.0 

Dl=2.*NMX+3. 

D2=Dl*(2.*NMX+5.) 

D3=D2*(2.*NMX+7.) 

D4=D3*(2.*NMX+9.) 

C     Using  Truncated  Taylor's  series  for  relative  upper  starting 

C      values  of  J(NMX)  and  J(NMX+1)  with  NMX  .GT.  ABS(Z*Z). 

C 

F1=1.-Z2/D1+Z2*Z2/(2.*D2)-Z2*Z2*Z2/(6.*D3) 
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F2=Z*(17D1-Z2/D2+Z2*Z2/(2.*D3)-Z2*Z2*Z2/(6.*D4)) 

C     Performing  downward  recurrence  with  highest  two  orders  in  CJ(N1). 
C 

M=NMX 
1 1  FSTO=Fl 

F1=(2.*M+1.)*F1/Z-F2 

F2=FSTO 

IF  (CABS(F1).LT.1.0E20)  GO  TO  1 

Fl=Fl*1.0E-20 

F2=F2*1.0E-20 

FSTO=FSTO*1.0E-20 
1  CONTINUE 

M=M-1 

IF  (M+2.GT.N1)  GO  TO  1 1 

CJ(N1)=F2 

CJ(N1-1)=F1 

N0=Nl-2 

C     Continuing  downward  recurrence  to  CJ(1) 
C 

DO22K=l,N0 
J=N1-K-1 

22  CJ(J)=(2.  *  J+ 1  .)*CJ(J+ 1  )/Z-CJ(J+2) 

C     Normalizing  entire  array  wrt  actual  J0(Z). 
C 

ALPHA=CSIN(Z)/CJ(1) 

D0  33K=1,N1 
33  CJ(K)=ALPHA*CJ(K) 

C     Computing  derivatives 
C 

DCJ(1)=CJ(1)/Z-CJ(2) 

D044K=2,N1 
44  DO(K)=a(K-l)-(K-l)*CJ(K)/Z 

RETURN 

END 

*********************************** 

SUBROUTINE  LPAD  (M,N1,C,P,DP) 
C 

C  This  routine  computes  Normalized  Legendre  polynomials  and 

C  derivatives. 
C 

C  For  m  >=  1;  P(n)  =  P(m,n+m-l)/sin(THETA) 
C  DP(n)  =  dP(m,n+m-l)  /  dTHETA    for  n=l,2,... 

C  Form  =  0;  P(n)  =  P(0,n)/sin(THETA) 
C  DP(n)  =  dP(0,n)  /  dTHETA        forn=l,2,... 

C 

REAL  P(*),DP(*) 
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KSP=N1-1 

Q=1.0-C*C 

S=SQRT(Q) 

IF(M.GT.O)  GO  TO  1 

P(l)=1.0 

P(2)=C 

DP(1)=0.0 

DP(2)=-S 

Q=S 

GO  TO  3 
1IF(M.EQ.1)G0T0  5 

P(l)=-((-S)**(M-l)) 

DP(1)=-M*C*(-S)**(M-1) 

GO  TO  7 
5P(1)=-1.0 

DP(1)=-M*C 
7D0  2K=1,M 

P(1)=P(1)*(2*K-1) 
2DP(1)=DP(1)*(2*K-1) 

P(2)=P(1)*C*(2*M+1) 

DP(2)=(C*DP(1)-Q*P(1))*(2*M+1) 
3D04K=2,KSP 

N=K+M-1 

AK=M-N-1 

P(K+ 1  )=((M+N)*P(K- 1  )-(2*N+ 1  )*C*P(K))/AK 
4  DP(K+ 1  )=((M+N)*  DP(K- 1  )+(2*N+ 1  )*(Q*P(K)-C*DP(K)))/AK 

RETURN 

END 

FUNCTION  ARCTAN2(Y,X) 
C 

C   Routine  returns  ATAN(Y/X)  for  all  cases 
C 

REAL  ARCTAN2,Y,X,PI 

PI=3. 14 15927 

IF(Y.NE.0.0)  GO  TO  1 1 

IF(X.EQ.0.0)  GO  TO  33 

IF(X.GT.0.0)  ARCTAN2=0.0 

IF(X.LT.0.0)  ARCTAN2=PI 

RETURN 
1 1  IF(X.NE.0.0)  GO  TO  22 

IF(Y.GT.0.0)  ARCTAN2=0.5*PI 

IF(Y.LT.0.0)  ARCTAN2=1.5*PI 

RETURN 
22  ARCTAN2=ATAN2(Y,X) 

RETURN 
33  ARCTAN2=0.0 

RETURN 

END 

p  ******************************************    ***************^******* 

C     PLOT3D  Plotting  subroutine 
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C     REVISED:  22  NOV  1989 

C 

C     This  subroutine  automatically  scales  and  plots  data  as  3-D  surfaces 

C 

C     All  subroutines  beginning  with  'Q' 

C  are  from  the  Grafmatic  plotting  library 

C 

SUBROUTINE  PLOT3D(Z,M,N,ASPECT,GRFLAB, 

+  TITLE  l,TrTLE2,HDRlyHDR2ALPHA,GAMMA) 

REALX(50,50),Y(50,50)Z(50^0)^MIN,XMAX,YMIN,YMAX^MAX2MIN 

REAL  XST,XFIN,YST,YFINZSTrZFINZFLOORZAVG 

REAL  Pffl,THETA^LPHA,GAMMA,STHETAJUl  ,M2,N2,R 

INTEGER  M.N 

CHARACTER*  1  YES  ,GR  AFPRNT 

CHARACTER*5ASTRING,GSTRING,TSTRINGJ>STRING>1IN,MAX,AVG 

CHARACTER*  12  name,label 

CHARACTER*64  TITLE  1,  TITLE2,  TITLE3,  TITLE4,  GRFLAB,  HDR1,  HDR2 

CHARACTER*64  HCOPY 

DIMENSION  IWORK1(640),IWORK2(640) 

C  SETTING  CONSTANTS 
C 

BELL=CHAR(7) 

FF  =  CHAR(12) 


C  Set  constant  for  degree  to  radian  conversion 
C 

PI  =  3.14159265359 

DTR=  PI/180. 

C  Calculate  #  of  points  on  array 
C 

NARRAY=PI*(  (M/2)**2) 

M2=FLOAT(M-1)/2.0 

N2=FLOAT(N-1)/2.0 

C    IWORK1  and  IWORK2  should  have  a  dimension  equal  to  the  number  of 

C    pixels  (or  points  on  a  pen  plotter)  in  the  x  direction. 

C 

NPTS=640 

C  NPTS  is  the  size  of  the  IWORK1  and  IWORK2  arrays 
C 

MODEP=16 

C  Plot  will  be  in  graphics  mode  16 
C 

C   XYSET  is  a  routine  called  to  set  up  the  initial  x,y  values  (see  below) 
C 

CALL  XYSET(X,Y,M,N) 
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c 

XST  =  X(1,1) 
YST  =  Y(1,1) 
ZST  =  Z(1,1) 

XFIN  =  XST 
YFIN  =  YST 
ZFIN  =  ZST 

ZMIN=1000 
ZMAX  =  -1000 
ZAVG  =  0.0 

PHI  =  45.0 
THETA  =  75.0 

DO10I=l,M 
D0  15J=1,N 
IF  (XG  J)  .LT.  XST)  XST  =  X(I  J) 
IF  (X(U)  .GT.  XFIN)  XFIN  =  X(IJ) 

IF  (Y(I J)  .LT.  YST)  YST  =  Y(IJ) 
IF  (Y(I  J)  .GT.  YFIN)  YFIN  =  Y(I  J) 

IF  (Z(I  J)  .LT.  ZST)  ZST  =  Z(I  J) 
IF  (Z(I  J)  .GT.  ZFIN)  ZFIN  =  Z(U) 

C  Check  if  point  is  on  the  antenna  array 

C 

I1=FLOAT(I)-1.0 

J1=FLOAT(J)-1.0 

R=  SQRT(  ((Il-M2)/M2)**2.0  +  ((Jl-N2)/N2)**2.0 ) 

IF(R.LT.1.0)THEN 

IF  (Z(IJ)  .LT.  ZMTN)  ZMIN  =  Z(IJ) 
IF  (Z(U)  .GT.  ZMAX)  ZMAX  =  Z(IJ) 
ZAVG  =  ZAVG  +  Z(IJ)/(FLOAT(NARRAY)) 
ENDIF 
15       CONTINUE 
10    CONTINUE 
XFTN  =  XFTN  +  .2 
YFTN  =  YFTN  +  .2 


C   The  following  parameters  set  the  axes  values,  etc.  (see  documentation 

C    for  Q3DXAX,  Q3DYAX,  Q3DZAX) 

C 

XBEG=XST 

YBEG=YST 

XEND=XFTN 

YEND=YFTN 

XSTEP  =  .2 
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YSTEP=.2 

IF  ((XFIN-XST).GE.  0.9)  XSTEP  =  .5 

IF  ((YFIN-YST).GE.  0.9)  YSTEP  =  .5 

ZDEF=ZFIN-ZST 

IF  (ZDIF.LE.  0.3)  THEN 

ZSTEP  =  0.05 

ZFIN  =  FLOAT(INT(10*ZFIN))/10+2*ZSTEP 

ZST  =  FLOAT(INT(10*ZST))/10-ZSTEP 

ZMINOR  =  1 
ELSEIF  (ZDIFXE.  1.0)  THEN 

ZSTEP  =  0.25 

ZFTN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 

ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 

ZMINOR  =  1 
ELSEIF  (ZDIF.LE.  2.0)  THEN 

ZSTEP  =  0.5 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 

ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 

ZMINOR  =  1 
ELSEIF  (ZDIFLE.  5.0)  THEN 

ZSTEP  =1.0 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 

ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 

ZMINOR  =  1 
ELSEIF  (ZDIF.LE.  10.0)  THEN 

ZSTEP  =  2.0 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 

ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 

ZMINOR  =  1 
ELSE 

ZSTEP  =30.0 

ZFIN  =  INT(ZFIN+ZSTEP) 

ZST  =  INT(ZST- 15.0) 

ZMINOR  =  1 
ENDIF 
ZFLOOR=0.0 


C   Clear  the  screen  and  enter  input  data 
C 

CALL  QSMODE(2) 

C  3D-stick  is  a  "wire-frame"  figure  with  hidden  lines  removed. 

C  3D-fill  is  similar,  except  that  the  interior  regions  of  each 

C     "panel"  making  up  the  surface  if  filled  in  with  color  "klrin" 

C    and  the  boundary  of  each  panel  is  drawn  with  color  "klredg". 

C  NOTE:  The  two  routines  use  different  hidden  line  removal  algorithm. 

C        Which  you  use  is  a  matter  of  taste.  However  the  equivalent 

C        of  Q3DFTL  is  not  available  on  the  pen  plotter  since  it  is 

C        based  on  an  overpainting  scheme  which  isn't  practical  on  a  plotter. 
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c 

C  The  angles  phi  and  theta  specify  the  view  direction  as  shown  in  the 

C        appropriate  figure  in  your  GRAFMATIC/PLOTMATIC  manual. 

C  The  parameter  icross  is  one  if  you  wish  cross-grid  lines  in  both  the 

C        x  and  y  directions  or  zero  for  grid  lines  only  parallel  to 

C        the  x  axis  (as  may  be  approrpiate  in  creating  a  spectral  plot). 

C 

C 

16  WRITE(*,*) '  Enter  0  for  3D-stick  or  1  for  3D-fill  option' 
C  REAEKV)  IFIL 

IFIL=0 

WRITE(*,*) 

WRITE(V) '  Default  view  angles  are:  phi=-45.,theta=70." 

IF(IFIL.EQ.1)THEN 

C  KLRIN  is  the  color  inside  the  panel 

C 

KLRIN=2 

C  KLREDG  is  the  panel  border  color 

C 

KLREDG=1 
ELSE 
C  ICROSS  is  1  for  cross-grids  in  both  directions 

C 

ICROSS=l 
ENDIF 

C   Continuation  point  to  loop  through  desired  angles  of  observation 
C 

17  CONTINUE 

CALL  QSMODE(MODEP) 

C    Modify  THETA  to  reflect  YFTN/ZFIN 
C 

STHETA  =  ATAN(  (YFTN/ZFTN)*TAN(THETA*DTR)  )/DTR 


C   Q3DROT  changes  the  view  perspective  as  described  in  the  documentation 
C 

CALL  Q3DROT(X,Y,Z)M,N,PHI,STHETA) 
C 

C ALL  Q3DWTND(XST,XFIN,YST, YFINZSTZFTN,XMIN,XMAX , YMIN, YMAX) 

IOPT=0 

IF  ( (XMAX-XMIN).NE.0.0)  THEN 
YOVERX  =  (YMAX-YMTN)  /  (XMAX-XMTN) 

ELSE 
YOVERX  =  999 

ENDIF 

XORG=0.0 

YORG=0.0 
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CALLQPLOT(100,550,30,270,XMIN,XMAX,YMIN,YMAX) 
1  XORG,YORG,IOPT,YOVERX,ASPECT) 


C  The  next  3  routines  create  the  x,  y,  and  z  axes 

C  If  Q3DFTL  is  used,  then  draw  the  axes  AFTER  the  3D  plot  is  drawn 

C  so  that  the  3D  plot  does  not  overwrite  the  x/y  axes. 

C 

IF(IFIL.EQ.l)  GO  TO  66 

C 

CALL  Q3DX  AX(XST^CFTN,XSTEP,9, 1 ,2,YBEG,YENDZFLOOR,  1 .0) 
CALL  Q3DY  AX(YST,YFTN,  YSTEP.9, 1 ,2,XBEG  .XEND.ZFLOOR,  1 .0) 
CALL  Q3DZAX(ZST,ZFIN,ZSTEP,9,-  1,2,XBEG,YBEG,  1 .0) 

66     CONTINUE 


C  Finally  draw  the  3D  plot,  using  the  appropriate  option. 

C 

IF(IFIL.EQ.l)  CALL  Q3DFTL(X,Y,M,N,KLRIN,KLREDG) 

IF(IFIL.NE.l)  CALL  Q3DSTk(X,Y,M,N,IWORKl,rWORK2,NPTS,ICROSS) 

C    Draw  the  axes  if  using  Q3DFIL  routine 

C 

IF(IFILJSE.l)  GO  TO  67 

CALL  Q3DX  AX(XST,XFIN,XSTEP,9, 1 ,2,YBEG  ,YENDZFLOOR,  1 .0) 
CALL  Q3DY  AX(YST,YFTN,  YSTEP.9, 1 ,2,XBEG  .XENDZFLOOR,  1 .0) 
CALL  Q3DZAX(ZST,ZF1N,ZSTEP,9,- 1 ,2,XBEG,YBEG,1 .0) 
67  CONTINUE 

Q  ********************* * ******************************************** 

C     WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  24 
CALL  QPTXT(NCHAR,TrTLEl,KOLOR(ICOL,IROW) 

Q  I,***************************************************************** 

C     WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  23 
CALL  QPTXT(NCHAR,HDRl,KOLOR,ICOL,IROW) 

C     WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  22 
CALL  QPTXT(NCHAR  ,HDR2,KOLOR,ICOL,IROW) 

C      WRITE  TEXT  ON  SCREEN 
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NCHAR  =  64 

KOLOR=10 

ICOL  =  10 

IROW  =  21 

CALL  QPTXT(NCHAR,GRFLAB,KOLOR,ICOLJROW) 

Q  ***********+************************************+*************+**,■ 

C     WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  15 
IROW  =  20 
CALL  QPTXT(NCHAR,TnLE2JCOLORJCOL,IROW) 

Q  ****************************************************************** 

C     MAGNITUDES  INTO  A  STRING 
C 

CALL  NUMSTR(ZMIN,MIN) 

CALL  NUMSTR(ZMAX,MAX) 

CALL  NUMSTR(ZAVG,AVG) 

TITLE3='min.  Z  =  '//NUN/A  max.  Z  =  '//MAX 

+    /A  average  Z  =  '//AVG 
C     WRITE  TEXT  ON  SCREEN 
C 

NCHAR  =  64 

KOLOR  =  10 

ICOL  =  17 

IROW  =  3 

CALL  QPTXT(NCHARTITLE3,KOLOR>ICOL,IROW) 

Q  ************** **************************************************** 

C     ANGLES  INTO  A  STRING 
C 

CALL  NUMSTR(ALPHAASTRING) 

CALL  NUMSTR(GAMMA.GSTRING) 

CALL  NUMSTR(THETATSTRING) 

CALL  NUMSTR(PHI,PSTRING) 

TITLE4=lnc.  A  =  7/ASTRING//,;Plane  G  =  7/GSTRING 

+    //;  View:  phi=,//PSTRING//'  theta='//TSTRING 
C     WRITE  TEXT  ON  SCREEN 
C 

NCHAR  =  64 

KOLOR  =  10 

ICOL  =  10 

mow  =  2 

CALL  QPTXT(NCHAR,TITLE4,KOLOR,ICOL,IROW) 

Q  ****************************************************************** 

C      Hardcopy  Query 
NCHAR  =  30 
KOLOR  =  12 
ICOL  =  1 

mow=i 

HCOPY='Hardcopy  — >  Enter  P  or  p' 

CALL  QPTXT(NCHAR,HCOPY,KOLOR,ICOL,IROW) 

CALLQCMOV(32,l) 
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CALL  QINKEY(IEXTEND,IKEY) 

NCHAR  =  40 

KOLOR  =  1 

ICOL  =  1 

IROW=l 

HCOPY= 

CALL  QPTXT(NCHAR^COPY,KOLOR,ICOL  JROW) 
C    If  P*  or  'p'  then  prtsc 
C 

IF(  IKEY.EQ.80  .OR.  IKEY.EQ.l  12  )  THEN 
CALL  QPSCRN 

ENDIF 

q      ****************************************************************** 

CALL  Q3DINV(X,YtZ,MJST) 

CALL  QSMODE(2) 

WRITE(*,*)ENTER  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE' 

RE  AD(5  ,*)Pffl,THETA 

IF  (  (PHI.NE.9.).AND.(THETA.NE.9.)  )  GOTO  17 
999     CONTINUE 
lOOFORMAT(Al) 
520  FORMAT(3X,I5) 
578  FORMAT(3X312.6) 

101  FORMAT(A) 
RETURN 
END 

***************************** 

Q     ************************************************************* 

c 

SUBROUTINE  XYSET(X,Y,M,N) 
C 

REAL  X(M,N),Y(M,N),IlrJl,N2,M2 
C 

C  define  x,y  values  over  the  m  by  n  grid 
C 

M2=FLOAT(M-1)/2.0 

N2=FLOAT(N-1)/2.0 

DOHlI=l,M 

D0  121J=1,N 

I1=FLOAT(I)-1.0 

J1=FLOAT(J)-1.0 

X(LJ)=  (Il-M2)/M2 

Y(IJ)=(J1-N2)/N2 
121        CONTINUE 
111        CONTINUE 

RETURN 

END 

p*  ********************** ********** 

C 

C  SUBROUTINE  NUMSTR 

C 
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C    ACCEPTS  NUMBER  AND  RETURNS  CHAR*5  WITH  SIGN  AND  THREE  DIGIT  STRING 
C 

SUBROUTINE  NUMSTR(ANUM.STR) 

REAL  ANUM.VNUM 

INTEGER  NUM,HUN,TEN,ONE,ORDER 

CHARACTER*5  STR 

ORDER  =  0 

C    Get  sign  information 
C 

IF  (ANUM.GE.0.0)  THEN 

STR(1:1)=V 
ELSE 

STROrlK-' 
ENDIF 

C    Strip  sign  values 
C 

VNUM=ABS(ANUM) 
C 

IF  (VNUM  LT.  1.0)  THEN 

ORDER=0 

VNUM= VNUM*  1 000.0 
ELSEIF  (VNUM  LT  .10.0)  THEN 

ORDER=l 

VNUM= VNUM*  100.0 
ELSEIF  (VNUM  LT.  100.0)  THEN 

ORDER=2 

VNUM= VNUM*  10.0 
ELSE 

ORDER=3 
ENDIF 

C   Use  integer  part  of  shifted  VNUM 
C 

NUM=VNUM 

IF  (NUM.GE.  1000)  THEN  NUM=999 

HUN=INT(NUM/100.) 

TEN=(NUM-100*HUN)/10 

ONE=(NUM-  100*HUN- 10*  TEN) 

C  PLACE  DECIMAL 
C 

IF  (ORDER  .EQ.  0)  THEN 

STR(2:2)=7 

STR(3:3)=CHAR(48+HUN) 

STR(4:4)=CHAR(48+TEN) 

STR(5:5)=CHAR(48+ONE) 

ELSEIF  (ORDER  .EQ.  1)  THEN 

STR(2:2)=CHAR(48+HUN) 

STR(3:3)='.' 
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STR(4:4)=CHAR(48+TEN) 
STR(5:5)=CHAR(48+ONE) 

ELSEIF  (ORDER  .EQ.  2)  THEN 
STR(2:2)=CHAR(48+HUN) 
STR(3:3)=CHAR(48+TEN) 
STR(4:4)='.' 
STR(5:5)=CHAR(48+ONE) 

ELSE 

STR(2:2)=CHAR(48+HUN) 
STR(3:3)=CHAR(48+TEN) 
STR(4:4)=CHAR(48+ONE) 
STR(5:5)=7 

ENDIF 

C 

RETURN 
END 
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Appendix  D 

SOFTWARE  SOURCES 

Microsoft  FORTRAN 

16011  NE  36th  Way 

Box  97017 

Redmond,  WA  98073-9717 

West  Coast  Consultants  CURVE  DIGITIZER 
4202  Genesee  Avenue,  Suite  309 
San  Diego,  CA  921 17 
(619)565-1266 

Microcompatibles,  Inc.  GRAFMATIC 
301  Prelude  Drive 
Silver  Spring,  MD  20901 
(301)593-5151 

Prof.  M.A.  Morgan 

Code  62Mw 

Naval  Postgraduate  School 

Monterey,  CA  93943 
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